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I. INTRODUCTION 


A. BACKGROUND 
1. A Carrier Navy 

The carrier is the nucleus of today's U.S. Navy, but its ability to project power 
is a byproduct of. the synergistic relationship it shares with its air wing. As Peter 
Garrison stated in CV: Carrier Aviation, "bv adopting submarine-launched ballistic 
missiles and ship-based tactical aircraft the navy has co-opted some of the tools of the 
army and air force and adapted them to a marine environment. Ships in themselves are 
no longer powerful, they take their power from the airplanes or long range missiles that 
thev carry. An adversary setting out to sink a ship does not send another ship against 
it, but instead sends an airplane or a missile” [Ref. 1]. (italics mine). The fleet has 
become dependent upon aircraft to take the battle over the horizon. and to extend the 
battle group’s area of influence. Without downplaying the importance of the carrier, it 
is just another large ship without its contingent of aircraft. - Thus, the carrier's real 
importance to the fleet rests on the ability of its personnel to support and maintain the 
air Wing at an optinial level of readiness. 

2. Aircraft Readiness 

Ideallv, the fleet would like each aircraft to be fully mission capable (FMC); 
that is, the material condition of the aircraft should satisfy minimum requirements to 
perform all of its missions. An aircraft's mission status (e.g., partially mission capable 
(PMC), not nussion capable (NMC), etc.) is determined using the Mission Essential 
Subsystems Matrices published in OPNAVINST 5442.4H. 

The ability to ensure a high percentage of FMC aircraft is constrained by 
several factors. First, the Navy has not been able to. provide shipboard intermediate 
maintenance activities (IMA) with sufficient numbers of qualified personnel, especially 
in critical technical billets. In part, this results from the allocation of manpower based 
on squadron size (1.e., the smaller squadrons have fewer personnel to maintain their 
aircraft). The problem associated with the manpower shortage 1s compounded bv a 
need for specialization in troubleshooting and repairing the different aircraft tvpes, 
models, and series deployed on the carrier. Second, the IMA’s capability to repair 


components as thev fail is constrained bv limitations of test equipment, maintenance 


facility size, availability of repair parts, and the total demand within the repair pipeline. 
Finally, those components which cannot be fixed at the intermediate level must be 
resupplied by outside sources. 

The air wing must be able to accurately identify and communicate its 
requirements for component replacements and resupply facilities should be capable of 
expeditiously responding to the those needs. During periods of high tempo flight 
operations, as would be experienced during wartime, the factors listed above become 
more critical. The ability to maintain high levels of readiness is hampered by an 
increase in the number of components in the repair pipeline, manpower and facility 
limitations, and a reduced or nonexistent ability ‘to resupply vital components. 
Potential problems may be partially rectified by improving the methodologies currently 
used to attain that “best” possible mix of replacement parts located on the aircraft 


Crier. [Ref 2: pp. 3.9] 


B. AVIATION CONSOLIDATED ALLOWANCES 
1. Definition | | 

The Aviation Consolidated Allowance List (AVCAL) is supposed to be a 
“oest’ mix inventory. An AVCAL is defined in OPNAVINST 4790.2C as a list of 
aeronautical material tailored to each individual ship to support assigned or embarked 
aircraft flight operations. This list of parts is prepared by the Navy’s Aviation Supply 
Office (ASO) under the direction of the aircraft type commander, and determines spare 
part allowances based on expected demands. 

Ideallv, it would be advantageous to stock at least one spare for each part 
that could possibly fail during deployment regardless of how unlikely that occurrence 
might be. Unfortunately, inventory sizes are ‘constrained bv both budget 
considerations and storage limitations. The latter is the result of a diverse collection of 
aircraft types, manufactured by different companies with each aircraft requiring its own 
specific replacement parts, coupled with the inherent aversion to building aircraft 
carriers even larger than thev already are simply to carrv more parts. Table 1 lists the 
normal complement of aircraft aboard a carrier, and is provided to illustrate the 
differences in aircraft types, missions, and spare part requirements. [Ref. 2: p.5]. 

2. Need for Improvement 
Determination of the proper AVCAL allowance for each component is not an 


easy task, and 1s definitely an area needing improvement. Á recent study conducted by 


TABLE 1 
NORMAL COMPLEMENT OF CARRIER AIRCRAFT 


AIRCRAFT MISSION SQUADRONS AIRCRAFT 
F-14A Fighter 2 24 
A-TE Light Attack 24 
A-6E Attack 10 


KA-6D Tanker 4 
S-3A ASW 10 
E-2C SURV 

EW 

ASW 





the Center for Naval Analyses reported that “statistics for combat aircraft operating 
during the period July 1982 through June 1983 show aircraft operating at less than 60 
percent of their anticipated wartime rate. These aircraft were not EMC due to supply 
between 18-24 percent of the time. These readiness rates would drop significantly if the 
aircraft flew at a wartime rate and AVCALs continued to be constructed as thev were 
in this period” (Ref. 3: p.1-1]. This statement strengthens the argument that current 
methods of AVCAL formulation are, at best, adequate for peacetime operations, but 
could result in disastrously low readiness rates when exposed to the additional 
requirements associated with wartime utilization. Two approaches to determining a 
better AVCAL would be the introduction of new methods of modelling requirements or 
the improvement of weak areas within the current model. 
3. AVCAL Models 

Two categories of models used to compute the AVCAL are readiness-based 
and demand-based inventory systems. Readiness-based systems determine spare 
allowances that achieve a predetermined level of readiness while minimizing the cost of 
the AVCAL. Examples of readiness-based models are Rand's DYNAMETRIC model 
and the Center for Naval Analyses MIME model [Ref. 3: p.2-S]. Demand-based 


(supply effectiveness) models determine component allowances based on expected 
demand and the capabilities of the repair facility and resupply pipeline. The ASO 
Manual model is an example of a demand-based method, and is the model currently 
used by the Navy to compute AVCAL spare allowances. 


C. ASO MANUAL MODEL 
The ASO Manual Model uses historic data to forecast demand rates and 
component repair times. The model determines the type and number of components 
included in the AVCAL allowance. For a detailed explanation of the rules governing 
component selection, the reader should refer to either Aviation Parts Allowance Policy 
by Peter Evanovich or A Retail Inventory Model for Naval Aviation Repairable [tens by 
Mark L. Mitchell. [Refs. 3,4]. 
l. Attrition Only Items 
Attrition only items are repairable components that cannot be repaired at the 
intermediate maintenance activity aboard the carrier. The expected demand for items 


of this type 1s computed using equation 1.1. 


NB’ = NB X (FH’/ FH) (eqn 1.1) 
where, . 
NB' = number of expected attritions for a specific flight program 
NB = number of attritions in the ASO data base 
FH' = specified flight program 
FH = total flight hours represented by the ASO data base 


2. Items that can be Repaired Locally 
The allowance for items that can be repaired locally is composed of two parts. 
An allowance for repairable components that are beyond the capabilities of the local 
maintenance activity to repair them (BCM), and an allowance for components that are 
queued at various places in the local repair cycle. 
a. BCM Allowance 
Units are considered BCM if the maintenance activity cannot repair the 
item because of a lack of tools, equipment, parts or technical expertise. Equation 1.1 1s 
used to forecast the number of BCM failures expecica to occur for a specific flight 


program. 


b. Repair Pipeline Allowance 
The second input to the local repairable allowance provides units to 
compensate for the number of failed components undergoing repair at the IMA. This 
allowance, called the Local Repair Cycle Allowance (LRCA), is set at a quantity which 
should ensure that demand does not exceed available supply at least ninety percent of 
the time. This safety level is dependent upon an accurate estimate of the number of 
components repaired by the shipboard intermediate maintenance activitv. Such a 

forecast is expressed in equation 1.2. 


NR'=NR * (FH'/FH) (eqn 1.2) 


Where, 


NR’ = number of expected repairs for a specific flight program 
NR number of repairs represented in the ASO data base 
and, 


.FH' and FH remain as defined above 


D. FAILURE RATE DISTRIBUTION 
Combining equations 1.1 and 1.2 provides the expected number of total failures, 
both repairables and attritables, for the specified flight program.! 


NR'+NB' = (NR* NB) X(FH'/FH) (eqn 1.4) 


Substituting NF’ for NR'+ NB' and NF for NRF NB, equation 1.4 can be expressed 
as: 


NF = NE SEM EH) (eqn 155 
Equation 1.6 results from rearranging terms on the right hand side of equation 1.5. 


NF' < (NE/FH) x FH' (eqn 1.6) 


l For “attrition only items” NR equals zero. 
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The total number of expected failures NF' equals the product of some specific 
future flight program (FH') and the quantity NF/FH. The latter should be easilv 
recognized as the maximum likelihood estimate (based on the historical data base) for 
the failure rate lambda. Thus, NF’ is the expected value of a Poisson distribution with 
the failure rate parameter A and interval of observation (0,FH']. The use of the Poisson 
distribution implies that the underlving distribution for the time between failures is 
exponential. Certain properties of the exponential distribution simplify the 
computational problem, but what price is paid for this simplicitv? The remainder of 
this paper will suggest alternative distributions for describing the time between 
component failures and estimating the expected number of failures for a given flight 


program. 


E. SYNOPSIS 

Chapter II describes the methods and assumptions used to extract component 
failure information from a data base containing six squadrons of F-14A aircraft flight 
and maintenance records. Chapters III and IV discuss current and possible methods 
of modeling aircraft failure rates, and techniques used to derive maximum likelihood 
estimates of model parameters. Flight hour dependent models are presented in Chapter 
III and include the traditional Exponential model and a Weibull model. Chapter IV 
discusses a sortie-dependent Geometric model, and the sortie-flight hour Mixed model 
proposed in Distinguishing the Effects on Failures of Changes in Sortie Rate and Sortie 
Length by Robert A. Levy [Ref. 5]. The four different models will be compared in 
Chapter V using both graphical and quantitative measures of fit. Finally, Chapter VI 
will discuss the advantages and disadvantages of each model, along with the 


conclusions and recommendations. 
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II. DATA BASE 


A. ORIGIN OF DATA 

The initial data base was obtained from the Navy Maintenance Support Office 
(NMSO), Mechanicsburg, Pennsylvania via the Center for Naval Analyses, Alexandria, 
Virginia. The data set included the historical flight and maintenance records for two 
operational deployments of F-14A aircraft. The most recent deployment included two 
squadrons, VF-114 and VF-213, which operated off the USS ENTERPRISE (CVN-65) 
during the period of March 1984 through December 1984. This deployment became the 
“model” set, with 1ts flight and maintenance records providing the parameter estimates. 
for the study. 

The second, earlier deployment consisted of four squadrons, VF-11, VF-31 and 
VF-14/VF-32, which operated off the USS KENNEDY (CV-67) and USS 
INDEPENDENCE (CV-62), respectively, during the period of November 1935 through ` 
March 1984. In this study these four squadrons were treated as the “validation” set. 

Using the parameter estimates from the model set, four specific failure rate 
distributions were used to calculate the expected number of component failures which 
should occur during a deplovment similar to the one flown by the four squadrons 
aboard' the USS KENNEDY and USS INDEPENDENCE. The actual ٣8٤٦ 
component failures observed during the deployment of these squadrons was then 
compared with the expected number generated using the parameters estimated from the 


“model” set deployment. 


B. MAINTENANCE DATA SYSTEM (MDS) 
|. Purpose f | 
The maintenance data system (MDS) was developed to improve the 
management capabilities of the Navy's Maintenance and Material Management (3- M) 
system. It was designed as a management information svstem which would provide the 
statistical data needed to analvze: 
e Equipment maintainability and reliability 


e Equipment configuration to include alterations and technical directive 
compliance 


e Equipment mission capability and utilization rates 


e Material usage 


e Maintenance and material procurement times 
e Weapon system and material costing 

To meet the data requirements listed above, NMSO divides the data collection 
effort into four specific areas: material reporting (MR), subsystem capability impact 
reporting (SCIR), maintenance data reporting (MDR), and utilization reporting. The 
material reporting system documents those materials used in support of the 
maintenance efforts. The subsystem capability impact report provides the higher level 
commanders with information pertaining to their subordinate command's ability to 
conduct their mission. The last two areas of data collection, maintenance data 
reporting and utilization reporting, provided the data necessary to determine the times 
between component failures. 

2. Maintenance Data Reporting (MDR) 

The documents which support the maintenance data reporting system include 
support action forms (SAF), metrology equipment recall (METER) cards, and the 
Visual Information Display System/ Maintenance Action Forms (VIDS; MAF). Support 
action forms document maintenance actions that did not require any corrective action, | 
such as, aircraft servicing, aircraft handling, and the preliminary look phase of an 
inspection. METER cards identify the maintenance time spent calibrating. and repairing 
the test and measuring svstems. The last form, VIDS'MAF (OPNAV 4790/60), 
documents inspections, technical directive compliance, and repair actions. It identifies 
the component failures by work unit code, time of failure (julian date and, when 
discovered code), aircraft experiencing the failure (bureau number), cause of failure 
(malfunction description code), and type of maintenance required (transaction code and 
action taken code). 

3. Utilization Reporting 

The flight records (utilization rates) were compiled from information recorded 
on the Aircraft Yellow Sheet (OPNAV 3760/2B). The Yellow Sheet is prepared by air 
crew personnel at the end of each mission and contains information identifying the 


flight according to aircraft flown, mission type, julian date, flight duration, etc. 


C. SELECTED WORK UNIT CODES 
A work unit code (WUC) is an alphanumeric code identifying an item on which 
work is being performed. Svstems are identified with two digit codes and components 


with five digit codes. Seven digit codes are used to further breakdown components into 


I» 


lower levels of subassembly. The components selected for this study are listed in 
Table 2. Each of the nine components chosen had at least thirtv failures during the 
deployment cycle. 


TABLE 2 
COMPONENTS WITH ASSOCIATED WORK UNIT CODES 


COMPONENT WUC 
CN1263/ASN92(V) INERTIAL MEASURING UNIT 0 
ECU74/A SIGNAL DATA CONVERTER 40 
MXU611/A FUEL TANK RELEASE MECHANISM 4622100 
CP1166/A AIR DATA COMPUTER 56X2500 
AN/ARCIS9/159(V)-1 UHF TRANSMITTER 632Z100 
T1224/AWG9 RADAR TRANSMITTER 74A 1500 
IP1185/AWG9 DETAIL DATA DISPLAY 7+40 
CP1248/ASW43 AIRCRAFT ROLL COMPUTER | 5772200 
1P1027/AVA12 ANALOG DISPLAY INDICATOR 6918100 





D. DATA REDUCTION 

The intent of the data reduction process is to eliminate those maintenance 
actions which do not identify a component failure (e.g.. aircraft servicing, 
troubleshooting, calibrations, or other routine maintenance actions). The maintenance 
codes not indicative of a component;item failure are found within the series of 
malfunction description, transaction, and action taken codes. [Ref. 6: pp. 17-19] 

1. Malfunction Description Codes 

The three digit numeric malfunction description codes describe, as accuratelv 

as possible, the probable cause of a system/component failure. Malfunction 
description codes used when no failure was detected are displayed in Table 3. 
Maintenance records which contained these codes were eliminated from the final data 
set. Conditional malfunction description codes identify component failures caused by 


factors other than natural wear-out. Examples of this tvpe of failure include: failed due 


to weather, foreign object damage (FOD), improper or faulty maintenance, etc. To 
reduce the complexity of parameter estimation, conditional failures were treated as if 
the failure had occurred through normal usage. The components studied had a low 
percentage of conditional failures. In fact, less than three percent of all observed 


failures were conditional. 


TABLE 3 
EXCLUDED MALFUNCTION DESCRIPTION CODES 


DESCRIPTION 


No defect 

No defect - removed and/or reinstalled to facilitate maintenance 
No defect - removed for modification 

No defect - removed for time change 

No defect - removed for scheduled maintenance 

No defect - removed for pool stock 

No defect - removed as Doo of a matched system 

No defect - removal directed bv higher authority 

No defect - removed during troublčshooting 





2. Transaction Codes 
Transaction codes are two digit numeric codes used to identify the tvpe of 
data reported bv the maintenance activitv. Transaction codes used in the maintenance 
reporting system are listed in Appendix A. Table 4 lists those transaction codes 
indicative of maintenance actions that required removal, installation, and/or repair of a 
defective component. Maintenance records associated with these codes were included in 
the final data set. All other maintenance records were deleted. 
5. Action Taken Codes 
Action taken codes are one digit alphanumeric codes which describe the 
maintenance action performed on the system, component or item. A complete list of 
action taken codes is shown in Appendix B. Action taken codes associated with the 
removal, installation, and/or repair of a defective component were included in the final 


data set and are shown in Table 3. 


TABLE 4 
TRANSACTION CODES USED IN STUDY 


DESCRIPTION 


EN work, including engines, involving non- 
e components/items documented as failed parts. 


repatrab 


Will be used when a repairable component is removed, 
excluding engine components, for processing at an I or D 
level maintenance activity (This code 1s used when 

only the removal must be documented and a replacement 
کا‎ not required). 


Removal and replacement of a defective or suspected, 
defective component from an end item, excluding engines 
at the O-level. Additionally, this code 1s used for the 
removal and replacement of a complete engine assembly 
for a, defect, suspected defect, or scheduled maintenance 
requirement. 


Will be used when a repairable engine component 
is removed for processing at an I or D level activity 
(This code is used when only the removal must be 
documented and a replacement 1s not required). 


Removal and replacement of a defective or suspected 
defective repairable component from an engine. 





E. FINAL DATA SET 
1. Problems 

Several problems were encountered during the data reduction process that 
require discussion. First, it was not always clear which flight evolution caused a 
component to fail. The cause for this confusion is the poor record keeping associated 
With the actual time of component failures. For example, the maintenance data 
reporting system does not record the actual time of component failure. Instead, the 
data is arranged according to the time the maintenance activity received the component 
for repair. In most, but not all cases, the component was received bv the maintenance 
activity immediately following the flight causing the component failure. 

Second, several records indicated multiple failures of a single component. For 
example, a single component may experience several unrelated malfunctions during a 
single flight evolution. To expedite the repair andor replacement of the failed 


component, the maintenance reporting system requires the documentation of each 


TABLE 5 
ACTION TAKEN CODES 


CODE DESCRIPTION 


1 BCM - maintenance activity not authorized to conduct 
repair 

2 BCM - Lack of equipment, tools, or facilities 
BCM - Lack of technical skills necessary to complete 
repair 

4 BCM - Lack of parts 


BCM - Fails check and test, and maintenance is allowed 
to conduct check and test only 


BCM - Lack of technical data 

BCM - Bevond the authorized repair depth 
BCM - Administrative | 
BCM - Condemned, repair not feasible 


ہے 3 مہ مہ C)‏ 


This code is used when a repairable item of material 
identified by a work unit code 1s repaired. 


"T1 


Failure of components; items undergoing test. 


"Y 


This code is entered when an item of material is 
removed and onlv the removal is to be accounted for. 


Q This code is entered when an item 1s installed and 
onlv the installation action 1s to be accounted for. 
R 


This code is entered when an item of material is 

removed due to suspected malfunction and the same or a 

like item is reinstalled. 
* BCM: Action taken codes labelled 1-9 are used for repairable 


items which have been administratively and/or technicallv screened 
and found to be bevond the capability of the maintenance activity. 


malfunction (1.e., multiple failures of a single component). Previous studies have 
treated this problem as independent single component failures with identical 
interoccurrence times. The alternative is to treat the multiple failure as a single 
incident. 

Third, the underlying assumption of the exponential distribution is that 
components fail as a function of flight hours flown. This would implv that a 


component is operated continuously, which may or may not be the case. 


tJ 
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Fourth, it was difficult to determine if a component was replaced immediately 
following the flight evolution causing the failure. Additionally, there was some concern 
pertaining to periods of flight operations during which a component might be missing. 

Fifth, the condition of the component (ie. new or used) at the start of 
deplovment was not alwavs known. A related problem involves the use of the same 
failure rate distribution for both repaired and new components. 

2. Simplifying Assumptions 

Before proceeding with any data analysis, it is necessary to make assumptions 
which can simplify or eliminate the problems discussed in the previous section. These 
assumptions are presented in the samie order as the problem areas they are intended to 
simplify. 

First, all failures are associated with the sortie during which they occurred for 
When Discovered Codes A, B, C, and D. For all other codes, it is assumed that the 
training and experience of the flight and ground crews ensure the discovery of the 
failure at the earliest- possible opportunity (1.e., the failure occurred during the flight 
prior to discoverv). A complete listing of When Discovered Codes is found in 
Appendix C. سے‎ 

Second, all multiple failures of a single component are treated as a single 
component failure. It is believed that multiple failures would affect the turnaround 
time required to repair the component and may cause delavs within the repair pipeline. 
However, this should not affect the rate at which the component fails. 

Third, it is assumed that components operate continuously, at least in a 
standby mode. If this is not the case, the sortie-dependent models should provide a 
Detter estimate Ol CECE allies: 

Fourth, the study assumes that components are replaced immediately and that 
gaps in coverage are not allowed. In fact, it is extremely unlikely that an aircraft 
would operate without any of the components listed in this thesis. 

Fifth, all components are considered new at the start of the deployment. For 
the Exponential and Geometric models this assumption is not necessary due to the 
“memoryless property” they exhibit. Additionally, the impact on the Weibull 
distribution should be minimal if the number of observed failures is large. 

Finallv, to reduce the complexitv of the problem, it is necessary to assume 
that repaired components exhibit the same failure characteristics as new components 


(1.e., have the same failure rate). 


3. Format 
a. Parameter Estimation 

To obtain the best possible estimate of the time between failures, the 
observations within the initial data base were sorted by aircraft bureau number, then 
Julian date and the event time. Failure times Were measured between consecutive 
component failures on a specific aircraft. These measurements were never made 
between component failures on two different aircraft. The information included in the 
final data set is the aggregate totals and includes a censoring indicator, number of 
sorties flown, flight time without failure, and flight time with failure. 

A censor indicator equal to zero indicates that a component failure was 
observed. The data entries associated with this type of observation include the number 
of flights flown inclusive of the event causing the failure, total flight time flown since 
last failure (does not include flight time of event causing failure), and the flight time of 
the event causing the failure. 

A censor indicator equal to one represents an event which is right hand 
censored. Right hand censoring refers to an observation of a component which does 
not fail prior to the termination of the deplovment (1.e., the data set documents a 
component s survival time, but provides no Information pertaining to the actual time of 
component failure). The data entries associated with this type of observation are the 
total number of flights flown and the total flight time since the last component failure. 

b. Model Comparison 

The intent of this thesis is to determine each model's ability to predict 
failures for the components listed, and compare those predictions with the actual 
component failures observed during the deplovment of the sguadrons aboard the 
۰۰۰۰۱۶۲۰۱۰ ظط‎ 1(۷ 1۱0 USS INDEPENDENCE. To accurately assess these capabilities, 
it was necessary to identifv component failures as they occurred during the deplovment. 
Ihus, the data set used to provide the graphical and quantitative comparison of the 
predicted vs actual component failures was different than the data set used to estimate 
the model parameters. The latter data set was sorted bv Julian date and event time 
only. Ihe accumulated number of sorties flown and flight hour totals were recorded 
with the observation of each failure, such that, comparisons could be made with the 


expected failures generated by the failure rate models. 


4. Summary Statistics 1 
Table 6 summarizes the two data sets and provides such information as the 


number of aircraft assigned, sorties flown, flight hour totals, component failures, etc. 


TABLE 6 
SUMMARY STATISTICS 


IE 


FLT HRS SONI 
DATA SET AIRCRAFT SORTIES FLIT HRS MEAN 0٦ 
MODEL 30 3489 7608.8 2.18 .98 
VALIDATION 44 3012 7030.4 27 .54 


f MODEL SET VALIDATION SET 
WUC FAILURES 0 


FAILURES 
734H 100 105 58 
46X 1600 34 49 
4622100 - 49 6 
56X2500 85 79 
6322100 49 31 
74A1500 107 75 
70 118 104 
5772200 . 34 33 
6918100 81 93 


MODEL SET = ٥+٥ 
VALIDATION SET < VF-ll, VF-l4, VF-3I, and VF-32 


III. FLIGHT HOUR MODELS 


A. OVERVIEW 

The two models discussed in this chapter describe component failures as a 
function of flight hours only. The first flight hour model to be examined is the 
exponential, which is the model currently used by the Navy to describe the distribution 
of the time between component failures. This model has some weaknesses inherited 
from the properties of the exponential distribution. The effect these problem areas 
might have ‘on the model's ability to accurately predict component failures will be 
discussed in the next section. The three models offered as alternatives should provide 
improvement to some or all of these problem areas, and the comparisons made in | 
Chapter V should shed new light on the magnitude of these deficiencies. 

The second flight hour model discussed is the Weibull. It is offered as an 
alternative method of describing the time between failures when flight hours are the 


only factors contributing to the component malfunction. 


B. EXPONENTIAL MODEL 
1. Weaknesses 

There are two major weaknesses associated with this model that cause some 
concern. First, the exponential is a flight hour model and as such it is insensitive to the 
effects of the sortie on the component’s time of failure. In Chapter V it will be shown 
that the expected number of component failures for the time interval (0,t] 1s 
proportional to t, so that the expected number of failures for a specific flight hour total 
is not affected by the average sortie length (i.e., the number of sorties flown). 

Second, the exponential distributions are “NO WEAR” distributions, which 
indicates that a component's probability of surviving a flight of duration t is the same 
for new components as it is for components which are "S" hours old (ie., the 
conditional survival probability is egual to the unconditional probability of survival). 
This is the "memoryless" property of the exponential distribution, and implies that the 
prior use of a comiponent has no effect on the time of its failure. 

The same concern is reflected in the exponential distribution's constant hazard 
rate r(t) (i.e., r(t)—A for all values of t). As such, the failure rate A would not be 


affected by the time on test t. It is believed that certain components exhibit hazard 


rates that are decreasing functions of time (t). In fact, many pieces of electronic 
equipment are routinely exposed to burn-in periods to make the component less 
susceptible to failure when put into operation. Conversely, it can be argued that other 
components have hazard rates that increase with time (1.e., a component is more likely 
to fail as 1t experiences some wear and tear). The exponential distribution is not 
capable of modelling these differences. 
2. Maximum Likelihood Results 
The maximum likelihood estimate, A', of the exponential failure rate 
parameter À (failures/hour) is found by iteratively solving equation 314 


Y T=Y t.exp(-At.) Tit ZO 0-0 am hr 


T OT bexpPCAu) — adc للا‎ 


(eqn sam) 


where, the left hand side of equation 3.1 equals the aggregate flight time for those 
sorties that did not experience a specific component failure, tj is the duration of the 
sortie causing the failure, n is the total number of records, censored and uncensored, 
within the final data set; and m is the number of uncensored observavions (failures). 

If the sortie lengths, : are equal to t for all j, then equation 3.1 can be 
simplified to the closed form expression given in equation 3.2. Thus, the average sortie 


length can be used to provide a “rough” estimate of the maximum likehhood estimator 


A: 


-In(FTWF/TFT) 
ee (eqn 3.2) 


where, FTWF (Flight Time Without Failure) equals the sum of the T;'s for i= l,...n, 
TFT equals the total flight time for the deplovment, and t equals the average sortie 
length. 

The maximum likelihood estimate for A was found using the FORTRAN 
program listed in Appendix E. The parameter estimates and 95 percent confidence 


limits for these estimates are displaved in Table 7. 


*For the derivation of this maximum likelihood equation, the reader is referred to 
Appendix D. 
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TABLE 7 
EXPONENTIAL MAXIMUM LIKELIHOOD RESULTS 


WUC M oh CL, 


734H100 0140 0018 (.0106,.0175) 
46X 1600 0048 0010 (.0029,.0067) 
4622100 0065 0011 (.0043,.0086) 
56X2500 0114 0021 (.0073,.0155) 
6327100 .0067 0012 (.0043,.0090) 
70 .0143 .0019 (.0105,.0181) 
74A5M00 .0159 .0024 ` (.0113,.0206) 
5772200 .0048 0011 (.0026,.0070) 
6918100 0108 - 0015 (.0079,.0137) 


K (I 2 PERCENT CONRIDENCE LIMITS: 





3. Quantile-Quantile Plots 

A graphical method often used to determine an empirical data's goodness of 
fit with respect to a theoretical distribution is the quantile:quantile plot. The technique 
requires the observations, Ap Ayres A, to be ordered from smallest to largest. 
Empirical quantiles are then computed for these order statistics based on a formula, 
such as Q.=(i-.5)/n, where Q, is the empirical quantile for the Xi) order statistic and n 
equals the total number of observations. These empirical quantiles are plotted against 
the theoretical quantiles to produce a graphical estimate of goodness of fit. 

Some adjustments to the above technique must be made to overcome the 
problems of censoring. The first adjustment relates to the problem associated with 
incomplete knowledge of the actual time of component failure. In Chapter II, it was 
mentioned that the failure occurred during the interval T. to T.+t.. To simplify this 
problem, it is assumed that the midpoint of that interval represents the actual time of 


failure. 
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The second adjustment that needs to be made accounts for those observations 
that are right hand censored. The Kaplan-Meier method for quantile estimation will 
account for those observations and will provide a more accurate estimate for the 
empirical quantiles. This method also requires that each observation, censored and 
uncensored, be ordered from smallest to largest. The empirical quantiles (Q.) are 


computed using eguation 3.3: 


Q; = 1-[(n+.5)/n] IT [(n-kK+.5)/(n-k+1.5)} for i in I (eqn 3°) 
k in I 
kK =1 
where, I is the set of all uncensored observations, and n is the total number of 
observations, censored and uncensored. [Ref. 7: p. 234] 

An exponential quantile-quantile plot was constructed for each of the nine 
components studied. The quantile-quantile plots for work unit codes 734H100 and 
56X 2500 are displaved in Figures 3.1 and 3.2, and the remaining quantile-quantile plots 
are presented in Appendix F. The y-axis of the plot represents the quantiles 0٤ 
and the x-axis represents the time between component failures. The asterisks identifv 
the uncensored observations, the X s identify: the locations of the censored 
observations, and the theoretical distribution evaluated using A” is defined by the solid 
line. The empirical data is fit well by the theoretical distribution if the data points are 
tightly packed around the solid line. Data points which deviate from the line can 
provide information about areas of weakness. ^ | 

Figure 3.1 illustrates the best exponential fit of all 116 11818 plots, 
but the plot still shows the lower to mid guantiles associated with the empirical data 
points overestimate the guantiles computed using the theoretical distribution. Figure 
3.2 provides a better representation of the other guantile-guantile plots. In every plot, 
the empirical quantiles overestimate the theoretical quantiles for quantiles less than .5, 
and underestimate the theoretical quantile for quantiles greater than .75. 

Figure 3.3 displavs the Weibull cumulative distribution function for & equal to 
5, 1, and 2. The reader should be aware that a Weibull distribution with shape 
parameter (4) equal to l is the exponential distribution. Figure 3.3 shows that for 
values of time less than to» the cumulative distribution for the curve represented by 
@=.5 is greater than the exponential curve (4=1), and for values of time greater than 


tg it is less than the cumulative distribution of the exponential curve. This is the 


31 


Exponential Quantile-Quantile Plot: WUC 734H100. 


Figure 3 | 
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identical situation displayed in each of the nine quantile-quantile plots. This would 
suggest that the distribution representing the time between component failures might 
be better modelled with a Weibull distribution with parameters, A> 0 and 0 Sa < 1. 


CUMULATIVE PROBABILITY 





a=WEIBULL SHAPE PARAMETER 


100 150 200 250 
TIME HRS 


Figure 3.3 Weibull Cumulative Distribution Function. 


C. WEIBULL MODEL 
|l. Properties 
The next flight hour model that will be discussed is the Weibull distribution, 
which can be a viable alternative to the exponential distribution since it is capable of 
modelling many different failure rates. The major diflerence between the Weibull and 
Exponential models is that the former is usually not “memoryless”, which indicates that 


the component’s previous usage can be exploited in the determination of the expected 


failure time. 
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The two parameters associated with the Weibull distribution are the scale 
parameter (A) and the shape parameter (4), with each shape parameter defining a 
specific failure rate distribution. As such, the Weibull family of distributions has the 
capability to model many different types of component wear. For example, 5 
greater than one describe a family of "WEAR OUT" distributions (increasing failure 
rate functions); an @ value equal to one describe the “NO WEAR” exponential; and 4 
values less than one describe a family of "WEAR IN” distributions (decreasing failure 
rate functions). Consequently, the Weibull distribution provides an additional 
modelling flexibility that should reduce the latter two areas of concern associated with 
the exponential distribution. 

The Weibull cumulative distribution function (Figure 3.3) and the exponential 
quantile-quantile plots indicate that a Weibull model with scale parameter (A> 0) and 
shape parameter (4 <1) might provide appropriate descriptors of the failure rate 
distribution for each of the components studied. As mentioned above, a shape 

parameter less than one is indicative of a “WEAR IN" family of distributions. 
Distributions of this type are characterized by a hazard rate function, r(t), which is 
decreasing in t, (Le., the component's failure rate decreases as t increases). Since mes 
majority of components studied are electronic components (AVIONICS), this 
assumption of a decreasing failure rate is probably an improvement over the constant 
failure rate of the exponential model. 

2. Maximum Likelihood Results 

Maximum likelihood solutions can be obtained by setting the partials equal to 
zero, and solving for 4’ and @’ using nonlinear optimization techniques sinular to the 
one described in Appendix G. The partial derivatives of the log-likelihood function 


with respect to the parameters A and @ are displayed below: 


7 a-l a T.U 
Koy E AY PY oque 
0۸ q کی‎ j i 
; ۳ | a 
81٤ ` y ACT; + ti] “logiA(T; + tp AT;) “log(AT Sp (ty) ` FAT) loe 
da < EN P | 


>The derivation of these equations is given in Appendix D. 


3 


ما 


where, 0,4 > 0; Tit 1,20; i7 1,2,...,m (the set of uncensored observations - failures); 
and j=m+1,m+2,...,n (the set of censored observations). 

The maximum likelihood estimates for A’ and q” were found using the 
FORTRAN program listed in Appendix G. The parameter estimates and 93 percent 


confidence limits for these estimates are displayed in Table 8. 


TABLE 5 
WEIBULL MAXIMUM LIKELIHOOD RESULTS 


WUC 7 CL, a CL, 


734H100 .0149 (.0084,.0216) .8079 (.8014,.8146) 
46X 1600 .0044 (.0006,.0082) .6664 (.6627..6703) 
4622100 .0066 (.0024,.0108) .7708 (.7666,.7750) 
56X 2500 0119 (.0063,.0176) .8170 (.8114,.8227) 
6322100 0067 (.0026,.0109) .7970 (.7928,.8011) 
7۸0 0165  '.  (.0084,.0246) .7066 (.6985,.7147) 
74A5M00 4 (.0104,.0224) .8844 (.8784,.8904) 
5772200 0045 (.0009,.0082) .7385 (.7349,.7422 
6918100 .0112 (.0058,.0167) 8142 (.8087,.8196) 


wai ٢ هده‎ 95 PERCENT CONPIDENCE LIMITS 





3. Quantile-Quantile Plots 

The parameter estimates listed in Table 8 were used to construct Kaplan- 
Meier quantile-quantile plots for each of the nine components studied. Two of these 
quantile-quantile plots are displayed in Figures 3.4 and 3.5, and the other seven are 
included in Appendix H. Both Figures graphically illustrate that the data is fit well bv 
the Weibull distribution, with the worst case fit represented by Figure 3.4. This plot 
shows that the Weibull distribution has a tendency to underestimate the empirical 
quantiles below ten percent, but fits the data extremely well for quantiles greater than 


ten percent. 


DS 


The results of the Weibull quantile-quantile plots are not conclusive, but offer 
credibility to the use of the Weibull model as an alternative method of modelling 


component failure rates when flight time is the only contributing factor. 
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IV. SORTIE DEPENDENT MODELS 


A. OVERVIEW 

The two models that will be discussed in this chapter differ from the models in 
Chapter III in that component failures can be affected by stresses incurred as a result 
of the sortie.* These stresses refer to events that occur during each flight, regardless of 
the sortie length, and would include events such as equipment initialization, takeoffs, 
landings, etc. The first sortie- dependent model to be discussed is the Geometric model. 
which determines a constant probability that a component fails on a given sortie. 

The second sortie-dependent model that will be briefly looked at is the Mixed 
model described in reference 5. This model could also be classified as a flight hour 
model since it determines the probability of a component failure based on both sortie 
and flight hour contributing factors. It is incorporated in this section because it also 
utilizes the geometric distribution to determine the probability that the first failure 


occurs-on the sth 


sortie. The major difference between these two models is that the 
Mixed model's probability of failure for a given sortie is not constant. In fact, it is a 
function of a constant failure probability caused by sortie-related stress, and a 


probability of failure which is dependent on the sortie length. 


B. GEOMETRIC MODEL 
1. Properties 

The basic premise of the Geometric model is that the component failure 
results due to sortie-related stress and is insensitive to differences in sortie length. This 
assumption can then be used to formulate the component’s survival or failure for a 
given sortie as a Bernoulli trial, where, the component fails with probability p or 
survives with probability 1-p. The geometric distribution is used to determine the 
probability of observing the first failure on the sth Bernoulli trial. 

The geometric distribution is in many ways similar to the exponential 
distribution and might be considered the exponential's discrete analog. Like its 
continuous counterpart, the geometric distribution is "memoryless" and has a constant 


failure rate, which would indicate that the Geometric model also ignores the history of 


*A sortie refers to one complete flight evolution, take-off to full-stop landing. 
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component usage. Any improvement offered by this model would be attributed to the 
effects of the sortie on the component's failure. 

It should be apparent that a low variability in sortie length (i.e., the sortie 
length is approximately constant) would imply that the exponential distribution is 
equivalent to a geometric distribution with p- l-exp(-At); where t equals the average 
sortie length. In this case, both models would provide similar estimates of component 
failures. The flight data used to generate the parameter estimates was acquired from 
aircraft flving during Flex-deck operations, and should have greater variability than 
would be expected from aircraft flying under normal cyclic operations." 

2. Maximum Likelihood Results 
The maximum likelihood expression for p', the estimate of the constant sortie 


failure probabilitv is given in eguation dane 


DES > s 2 ert (eqn 4.1) 


where, m represents the number of uncensored observations (sorties with component 
failures), and the sunimation of the sis for 15 l,2,...,n represents (hema! number of 
sorties flown. Thus, the maximum likelihood estimator p' equals the number of flights 
experiencing a component failure divided bv the total number of sorties flown. This is 
the same maximum likelihood estimate which would have resulted without censoring. 
The maximum likelihood estimates p' for p and the associated 95 percent confidence 


intervals are listed in Table 9. 


C. MINEDSMODEE 
1. Properties 
This section will provide a brief description of the Mixed model proposed bv 
Levy [Ref. 5]. This particular model is structured such that the sortie and the sortie 
length (flight hours) can combine to cause a component failure. It is obvious that 
factors contributing to one component's failure may or mav not be detrimental to 
another components operating life. Most likely, "the life of electronic components is 


probably affected more by the number of power surges to which thev are subjected 


"During Flex-deck operations, aircraft sortie length is not constrained to a one to 
two hour launch and recovery cycle. 


The derivation of this equation 1s given in Appendix D. 
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than operating hours; landing gear life should be influenced primarily by takeoffs and 
landings; and engine life should be affected by rotational speed changes, which are 
accompanied by heat transients and pressure changes” [Ref. 8: p. 2]. The Mixed model 
appears to have an increased flexibilty to model some of these sortie and flight hour 
factors. 

As in the previous section, a constant probability of failure is employed to 
model the effects of sortie induced stress. In addition, the Mixed model incorporates a 
failure probabilitv for those components whose lifetime is affected by the flight hours 
flown. Combining the constant sortie failure probability with the flight hour probability 
results in the expression for the probability of failure p; for the ¡th sortie, which is 
expressed as the complement of a component surviving the stresses of the sortie and 


the flight hours flown: 
Dee ne ED 0 CB <1, > 0, t 20 (eqn 4.2) 


where, the fight hour failure distribution is exponential and pg equals the constant 
sortie failure probabilitv. - 
2. Maximum Likelihood Results 
The partial derivative of the log-likelihood function with respect to the 


parameters À and p. are listed below:”: 


OK (1-pg)t;exp(- t) 


=. | xl. SET (eqn 4.3) 
Oh = 1-exp(-AtXI-po) 1 - J 

OK _ y exp(-At;) s;-l e S a 
Op, 1 l-exp(-At;)X1-p.) Epo j E l 


gere, OSp, Sl; À > 0: Tit; T; 1ل‎ cU msthe sef ol uncensored 
observations - failures); and }=m+1,m+2,....n (the set of censored observations). 
Maximum likelihood solutions for 4’ and p, were obtained using the 
FORTRAN program listed in Appendix G. For each of the nine components studied, 
the log-likelihood function was maximized when i’ equaled zero, and Po’ equaled p, 


the maximum likelihood estimate for the constant failure rate p of the Geometric 


The derivation of these equations is given in Appendix D. 


4] 


model. The results for A” are included in Table 9. The reader should note that the 
values for p' and pg' are identical, and the failure for À is zero for every component. 
Since both models provide the same maximum likelihood results, they will be jointly 


referred to as the Geometric model when compared with the other two models. 
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V. MODEL COMPARISONS 


A. PURPOSE 

The purpose of this chapter is to compare the Exponential model with the 
Weibull and Geometric models to determine if the distribution currently used provides 
the most accurate description of the failure process. As a method of comparison, the 
Chi-square goodness-of-fit test was not considered for two reasons. First, the failure 
rate information provided by the censored observations would have been lost because 
of the impossibility of determining the exact time of component failure. Second, 1t was 
considered more important to examine each model's capability to predict the number 
of component failures resulting from a specified number of sorties and flight hours 
flown. The next section describes the methods of estimating failures for each model. 
These estimates are then compared with the actual number of failures to determine 


each model’s goodness-of-fit. 


B. COMPONENT FAILURE ESTIMATION 
1. Geometric Model 
The determination of the expected number of component failures -associated 
with the Geometric model is straightforward. The parameter p estimated using the 
Geometric model defines the probability that a component would fail during a given 
sortie. Furthermore, if it is assumed that the outcome of a sortie is independent of the 
outcome of the other N-1 sorties, and p is constant for each sortie, then the probabilitv 
law defining the number of component failures for the deployment is the sum of N 
independent Bernoulli trials. This random sum has a Binomial distribution with an 
expected value equal to N Xp. Thus, using the Geometric model the expected number 
of component failures is calculated as the product of the number of sorties flown (N) 
and the probability that the component fails during the sortie. 
2. Exponential Model 
To compute the expected number of component failures using the Exponential 
model, each failure can be thought of as an occurrence within a Renewal process; a 
stochastic process which counts the number of occurrences within an interval (0,t]. It is 
assumed that the times between successive occurrences are independently and 
identically distributed non-negative random variables. As such. it can be shown that 


the expected number of occurrences within an interval (0,t] is equal to: 
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M(t) =>, Pv k= o m 
k 


where, M(t) equals the expected number of occurrences (failures) within the interval 
(0,t], and F,(t) is the k-fold convolution of the cumulative distribution function. The 


derivative of M(t) with respect to t defines the renewal rate of the process as: 


M(y=m(9)=) f (0)  k=1,2,3,... 
k 


Since the k-fold convolution of the exponential probability density function is known 
to be Gamma, m(t) is egual to the infinite sum of Ganima probability density functions 


with scale parameter A and shape parameter k. 
4 Kk,k-l,-Ax 
pe NS 
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m(t) 


In tLe equation above, DI (AE TA) for k=1,2,3,... 1s easily recognized as the 


Taylor's series expansion of eh 


Thus, the renewal rate for the counting process 
.defined by exponential interoccurrence times is A. The expected number of failures 


M(t) is found by integrating m(t) from 0 to t. 


t > 
M(t)= {A ds=At 
0 


Therefore, the expected number of component failures is equal to the product of the 
failure rate À and the time on test t. This particular Renewal process is referred to as a 
Poisson process. 
3. Weibull Model 

Estimating the number of expected failures associated with the Weibull model 
Is not as easy as the other methods discussed above. The observation of component 
failures is still a Renewal process, but the interoccurrence times are distributed 
according to a more complex distribution function. Finding the k-fold convolution of 
either the cumulative distribution function or the probability density function is an 


extremely difficult if not impossible task. For processes defined by continuous 
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interevent times with finite means and variances, the Elementary Renewal Theorem 
Zu? 


approaches infinity [Ref. 9: p. 290]. Unfortunately, this theorem does not help estimate 


states that M(t) is approximately Normal with mean- t/H and variance = t6 as t 
failures when t is small. As an alternate method to construct M(t) when t is small, a 
simple computer simuiation is used. 

The inputs required for the simulation included the number of aircraft 
assigned to the carrier air wings, the Weibull shape and scale parameters, and an 
approximate flight hour distribution. This information was used to repeatedly simulate 
the air wing operations during the deployment. 

The aircraft were chosen at random to fly a fictitious mission of some sortie 
length t. For the modelling set, there were 30 aircraft assigned. The model assumed 
that each aircraft had an equally likely (1/30 percent) probability of being chosen to fly 
a mission. Realistically, aircraft are not chosen at random to fly missions, and most 
likely those aircraft recently flown are better prepared for future flights. A weighting 
scheme could have been used to increase the probabilities of those aircraft flown 
recently, and decrease the probability for those aircraft that were undergoing siniulated 
maintenance activities. | | 

The duration of the mission was determined by a random variable generated 
from a Normal distribution with the mean and variance estimated by the Model data 
set (see Table 6). Figure 5.1 depicts a histogram of the sortie lengths for the Model set 
with an overlay of the appropriate Normal distribution. Between the values of Q to 6 
hours, the Normal distribution fits the histogram adeguately enough to justify its use 
as the sortie length's parent distribution. Values of randomly-generated sortie lengths 
less than .5 hours or greater than 6 hours were truncated to .5 hours and 6 hours, 
respectively. 

After the selection of aircraft, the conditional probability of failure was 
computed for a flight of sortie length t. The only information required for this 
computation was the component's Weibull shape and scale parameters. and the 
aircraft's flight history since the last component failure. The conditional probabilitv of 
failure could then be calculated as described in Appendix D. The value for the 
conditional probability of failure was then compared against a Uniform (0,1) random 
variable, and if the uniformly-generated outcome was less than or equal to the 


conditional probabilitv of failure, a component failure was simulated. 
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Figure 5.1 Sortie Length Histogram vs Normal Distribution. 


The number of component failures was recorded: at 500 hour intervals until 
the completion of the deployment. The deployment was simulated 100 times, and the 
average number of component failures was computed for each 500 hour interval. These 
averages were the approximations to the expected number of component failures M(t) 
at 500,1000,1500,... hours étc. 


C. COMPARISONS 
l. Methodology 
The predicted failures were compared with the actual failures as they occurred 
in each ol the two deployments (USS KENNEDY and USS ٢٦٢٦٢٢٣٠٠٦ 
this point, it is important to clarify the speeific calculations of the expected failures. As 
described earlier, the accumulated flight time and total sorties flown were recorded with 
the observation of each failure. Therefore, the exponential and geometric expected 


failures could be computed using the formulas listed im this chapter. Unfortunately, the 
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expected number of Weibull component failures was not always defined at the 
accumulated flight time of the actual failure. To provide such an estimate, it was 
necessary to interpolate? between the Weibull expected failures associated with the 
simulated failure times above and below the actual time of component failure. The 
difference between the expected and actual failures was computed at the time of each 
component failure. These differences were then summed and divided by the total 
number of component failures to provide an estimate of the average prediction error. 
An average difference close to zero, with small standard deviation, indicates that the 
model is a good predictor of actual failures. A negative average shows that the model 
has a tendency ‘to underestimate the number of component failures, and a positive 
average is indicative of overestimation. Tables 10 and ll provide the average 
differences and corresponding standard deviations for the “model” and “validation” data 
sets respectively. 

The model providing the best fit for a specific component was chosen based 
on the average difference between the predicted failures and the actual failures. Table 
12 provides a summary listing of the optimal distributions for each component in both . 
data sets. The reader is cautioned that the optimal was chosen based on an average 
difference closest to Zero, and in some cases the expected number of failures was 
adequately described by other distributions as well as the one listed as optimal. If two 
or more models provided similar estimates, the one with the smallest standard deviation 
was chosen as optimal. The actual failures were plotted against the predicted failures to 
further illustrate a distribution's ability or disability to estimate failures. These graphs 
are included in Appendix ٠ 

2. Results. 

Before proceeding with the comparisons, a caveat should be placed on the 
results. The abilitv of the models to predict the failures observed during the deplovment 
of the “validation” data set is contingent on the accuracy of the model parameter 
estimates. For a valid comparison, it is necessary that both deplovments come from the 
same population. That 1s to say, the flight operations were conducted under similar 
conditions. The reader should be aware that the “model” set aircraft and the 
“validation” set aircraft operated during different seasons and in different geographic 


locations. The parameters estimated from the “model” set may or may not provide the 


۱ > | کس 

SÅ quadratic interpolation of the form aFT" £ bFT £c< EF was used. Where. FT 
egualed accumulated Hight time of failure, EF egualed expected failures, and a,b,c were 
parameters estimated using the known Weibull times and expected failures. 
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TABLE 10 
MODEL SET AVERAGE DIFFERENES 


Exponential Geometric Weibull 

WUC :734H100 

Average: FO |د‎ 51 3420 

St.Dev.: 2 6 5.55 27883 

WUC :46X1600 

Average: 0.94 0250 3.65 

St.Dev.: 91 T 1398 
| WUC:4622100 

Average: - 2 OM RS 0.56 

SE Dent: 2.54 1. 2.52 

WUC:56X2500 

Average: . 21 7 07 1.43 

S وه‎ 1.94 1.16 1.79 

WUC:6322100 

Average: - -0.04 0.08 2.74 

ST DEV, SS 2 2 2:29 

WUC:74A1500 

Average: ٤ و‎ -0.03 578 

St.Dev.,: 4.77 3.92 5909 

WUC:74A5MOO 

Average: 220.53 14:05 0:08 

St.Dev.: 3.04 4.45 4.14 

WUC:5772200 

Average: -1.49 -1.84 1295 
OL DEN s 2:13 1:29 2.44 

WUC :6918100 

Average: -4.36 =2 97 Sq 

St.Dev.: 3125 2.53 3308 


true estimate of the population parameters required to accurately predict failures within 
the “validation” data set. 

Within the “model” set, all of the models predict failures adequately. This 1s 
not unexpected because the parameters defining the distribution Were estimated using 


the flight data from the “model” set. Table 11 shows that this is not true for the 
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TABLE 11 
VALIDATION SET AVERAGE DIFFERENCES 


Exponential Geometric Weibull 
WUC :734H100 
Average: 20.76 16,44 26.22 
St.Dev.: 11.01 8585 12201 
WUC :46X1600 
Average: aS 11 -8.04 -0.37 
5۲.۰067. : 4.04 5.352 4.17 
WUC:4622100 
Average: 16.20 14.57 20590 
SE Deva 15433 J 0.95 16:56 
WUC 0 
Average: -4.23 -7.41 0.48 
SC: Dev.: 4.01 2 5 3.49 
WUC :6322100 
Average: 25588 -8.43 2) 5 
SENA 3.14 4.45 2.42 ; 
WUC:74A1500 
Average: 5582 2.00 DEM 
5۲.067 ٤ 6810 4.26 8.06 
WUC :74A5M00 | 
Average: 8.67 2183 85 
5 ۰027.۰ 3.41 3 75 3.74 
WUC:5772200 
Average: -4.61 -6.33 27 
St.Dev.: pus 2.40 202 
WUC:6918100 
Average: -2 01 5-01 255] 
St.Dev.: 7:12 8.48 ILT 


deployment defined by the “validation” data set. In most cases, a single optimal model 
can be identified. 

Table 10 illustrates an extremely minor difference between the exponential and 
geometric estimates. In fact, it could be argued the geometric and the exponential were 


interchangeable as predictors of failures for the “model” set deployment. For the 
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TABLE 12 
OPTIMAL DISTRIBUTIONS 


WUC MODEL SEM VALIDATION SET 
734H 100 Exponential Geometric 

46X 1600 Geometric Weibull 

4622100 Weibull Geometric 


56X 2500 Geometric Weibull 
6322100 Exponential Weibull 

74A 1500 Geometric Geometric 
7۹۸۹0 Geometric Geometric 
57172200 Geometric Weibull 
0913100 Weibull Exponential 





Inertial Measuring Unit (WUC:754H100), the optimal distribution is listed as 
exponential for the “model” set and geometric for the "validation" set, but it could have 
been geometric in both cases. During the deployment represented by the “validation” 
data set the average sortie length increased by approximately 8.5 percent. As 
mentioned earlier, the exponential distribution is insensitive to changes in the average 
sortie length and would provide an estimate based on flight hours only. On the other 
hand, if the geometric distribution was correct the longer sortie lengths would imply 
that a specific flight hour total was accomplished with fewer sorties flown. It would 
provide a smaller estimate than would the Exponential model. For Work Unit Code 
734H100, the Exponential model was unable to account for the longer average sortie 
length (in the “validation” set) and overestimated the actual number of failures. 

The Geometric model appears to provide the best fit for components that are 
not necessarily operated continuously. Examples of components of this type would 
include the AWG-9 Radar Transmitter (WUC:74A1500) and Detail Data Display 
(WUC:74A5M00). A possible exception to this rule would be the Inertial Measuring 
Unit (WUC:734H 100). Despite its continuous use, the Inertial Measuring Unit's gvros 


and accelerometers are extremely sensitive to the movements and impacts associated 
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with carrier takeoffs and landings, and it is extremely likely that component failures are 
caused by sortie-related stress vice continuous flight operations. 

For continuous use components within the "validation" data set, the Weibull 
distribution provided the best estimate of component failures. The exponential 
distribution was the optimal predictor for only one component (WUC:6918100), and in 
that case, the Weibull distribution was nearly as accurate. From a readiness 
standpoint, the Weibull model might be considered a better choice because of its 
tendency to overestimate the actual number of failures. Amazingly, the exponential 
distribution predicted failures best for less than 10 percent of the components within 
the “validation” data set. In support of the current model, it should be noted that the 
exponential’s estimate of failures was more than adequate for the “model” set. 

As mentioned previously, the Weibull distribution with shape parameter equal 
to one is exponential. Figure 5.2 demonstrates the capability of the Weibull simulation 
to provide estimates similar to those generated for the Exponential model. This figure 
shows the Weibull and exponential predictions for the Detail Data Display 
(WUC:74A5M00) with estimated Weibull shape parameter equal to .91. ۱ 

It should not be surprising that in some cases, none of the models will be 
adeguate estimators of the actual failures. This type of phenomenon was observed for 
the Fuel Tank Release Mechanism (WUC:4622100). For two successive deplovments, 
total failures went from 49 to 6 for a similar number of flight hours and sorties flown. 
Based on the available data. it is hard to determine if 49 was too high or 6 was too 
low. This type of outcome emphasizes the difficulty in predicting the actual failures for 
an upcoming deployment based on parameters estimated from one deplovment rather 


than the entire historical data base. 
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VI. SUMMARY 


A. ADVANTAGES AND DISADVANTAGES 

This section provides a brief summary of the primary advantages and 
disadvantages associated with each of the four models discussed in Chapters III and 
IV. The major advantage of the Exponential and Geometric models is that the number 
of expected failures is easily calculated using closed form expressions. Both models 
express failures as a linear function of only one factor, flight hours (Exponential) or 
sorties flown (Geometric). There is a tendency then to concentrate on aggregate flight 
totals while ignoring the specific flight history of each aircraft. As a result, neither 
model is able to adjust its failure estimates to compensate for differences in flight 
programs. i 

The Mixed model, on the other hand, attempts to account for failures resulting 
from sortie stress and/or continuous use. Unlike the previous models, the expected 
number of component failures is 'a function of a specific flight program. The 
probability of failure associated with the Mixed model is a function of the sortie 
duration. As such, this model should be more sensitive to differences in flight: 
programs and would provide a more accurate estimate of failures than the Exponential 
and Geometric models. Use of the Exponential distribution to model continuous-use 
components has weakened the forecasting capabilities of the Mixed model. Lack of 
consideration for flight hours flown since the last component failure oversimplifies the 
problem. As stated previously, the small variability in sortie length combined with the 
“memiorvless” property of the Exponential distribution would implv that the probability 
of flight hour induced failure is approximately constant for all sorties flown. Since the 
Stress related failure probability is also constant, this model would not necessarily 
predict any better than the Geometric model. This tvpe of phenomenon was observed 
in this study. 

Ihe Weibull's strongest feature is its abilitv to model failure rates for components 
which have increasing, decreasing, or constant failure rates. Since it does not possess 
the memoryless property of the Exponential; Geometric models, the expected number of 
failures is dependent upon previous flight histories. Like the Mixed model, this should 


provide an improved estimate of component failures. One disadvantage is a lack of 


simple mathematical formulas with which to calculate the expected failures. However, 


in today s world of computers this should not be a major problem. 


B. CONCLUSIONS AND RECOMMENDATIONS 

It is extremely difficult to estimate the number of component failures which 
occur during a deployment based solely on the data provided by another deployment. 
Such an approach would not account for differences in flight and maintenance 
personnel, weather conditions, geographical location, mission type, aircraft's 
maintenance history, etc. For example, the number of actual failures observed for the 
Fuel Tank Release Mechanism went from forty-nine to six for a simular number of 
flight hours and sorties flown. To minimize the effects caused by such variations in 
deplovments it is recommended that these models be evaluated using an enlarged data 
base. 

Surprisingly, in this study, the model currently used by the Navy (Exponential) 


did not provide the best estimate of component failures. The Geometric provided a 


better fit for components which were not operated continuously, and the Weibull. 


performed better when the components were operated continuously. In fact, overall, 
the Exponential was the least effective model for these nine components. 

The results obtained. using the Weibull simulation were encouraging. The model 
produced satisfactory estimates with a minimal number of assumptions, data input, and 
software coding. While the results obtained using this model were not overwhelming, 
an improvement in prediction capability was observed. 

Two possible directions for improving the current method of estimating demand 
(1.e., component failures) are the construction of a universal model and the use of high 
resolution simulations. The first alternative suggests the use of a model which can 
describe components with different failure rate functions. and incorporates a 
methodology utilizing sortie-induced stress and flight hours flown to compute an 
estimate of component failures. The model currently used does not have this 
capability. 

A slight modification to the Mixed model described in Chapter [V nught provide 
the modelling flexibility desired. A failure rate distribution, such as the Weibull. could 
be used to describe continuous-use failures instead of the Exponential. Such a change 
would ensure that the component's failure probability could be influenced bv prior use. 


That is to say, the probability of failure for a given sortie would now be a function of 


the sortie length and the aggregate flight time since the last component failure. The 
strengths of the Weibull distribution combined with the constant probability of failure 
(sortie) associated with the Mixed model should provide an estimate of demand which 
is more sensitive to variations within and between flight programs. 

The price paid for this improvement in forecasting ability is an increase in the 
complexity of the expected demand computations. With the abundance of computers 
within the fleet this should not be a problem. Perhaps, ten to twenty years ago it was 
necessary to rely on simple models for prediction and the crude estimates they 
provided. This is not the case today; computer simulations can be used when simple 
mathematical formulas are not available. As described in Chapter V, a low resolution 
simulation was used to forecast the component failures associated with the Weibull 
renewal process. The advantage of such simulations is the abilitv to include other 
factors which may affect a component's failure rate. Examples of these factors would 
include number and type of landings, mission type, maintenance programs, weather 
conditions, etc. It is safe to say, that these more complicated, high resolution models 
could include more information about the factors which cause components to fail, and 


should ultimately result in better forecasts of demand. 
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APPENDIX A 
TRANSACTION CODES 


CODE DESCRIPTION 


Used to report an inventory gain. 


Used to report a change in the readiness reportable 
status of an equipment. 


Used to report an equipment loss. 


a. On-equipment work not involving a removal or 
replacement of a defective component or item. 

b. On supproting engine documents not having a removal 
of a defective or suspected defective component or item 
when the engine 1s not specifically identified to a 


‘ specific aircraft. 


c. Used at the O or I level maintenance activities when 
closing out a maintenance acuon. 


On-equipment work, including engines, involving non- 
repairable components:items documented as failed parts. 
Removal of a nondefective component, item, excluding 
cannibalization, from-an engine to be processed at the 
O-level maintenance activity. | 


Installation of a nondefective component item. 
eXcluding cannibalization, on an engine to be processed 
at an O-level maintenance activity. 


Removal of a nondefective component item, excluding 
engine components items and a cannibalization to be 
processed at an O-level maintenance activity. 


Installation of a nondefective component/item |. 
excluding engine component items and cannibalization. 


a. The removal and replacement of nondefective 

engines and component, items to accomplish a 
cannibalization action. 

b. The removal and replacement of those consumeable 
components; items subject to a scheduled removal interval 
or items of supply significance. i ۱ 

c. The removal and replacement of a nondefective engine 
component for cannibalizauon at the I-level only. 


a. The removal and replacement of nondefective 

engine component to accomplish a cannibalization at the 
O-level onlv. 

b. The removal and replacement of those consumeable 
engine components su ject to a scheduled removal interval 
Or items of supply signilicance. È 


Will be used when a repairable component is removed, 
excluding engine components, for processing at an | or D 
level maintenance activitv ( This code is used when 

only the removal must be documented and a replacement 
is not regutred). 


Removal and replacement of a defective or suspected | 
defective component from an end item, excluding engines 
at the O-level — f | 

Additionally, this code is used for the removal and 
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replacement of a complete engine assembly for a defect, 
suspected defect, or scheduled maintenance requirement. 


Will be used when a repairable engine component 
is removed for processing at an | or D level activity 
(This code is used when onlv the removal must be 
documented and a replacement 1s not required). 


Removal and II of a defective or suspected 
defective repairable component from an engine. 


Used to document components processed through the 
I-level maintenance activity for check, test, and 
service. 


Work performed on a removed repairable component with 
no failed parts or awaiting parts documented in the 
Failed/Required Material blocks of the VIDS; MAF. 


Work performed, on a removed repairable component with 

failed parts, awaiting parts, or cannibalization 

a documented in the Failed, Required Material 
ocks. 


Close out for man hours or awaiting parts at the 
I-level maintenance activity. 


. Technical directive compliance with no part number 


change. . 


Technical directive compliance with a part number 
change. 


Will be useed to report SCIR data bv the reporting 
custodian when transient maintenance is performed by 
other than the reporting custodian. 
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APPENDIX B 
ACTION TAKEN CODES 


CODE DESCRIPTION 


BCM - maintenance activity not authorized to conduct 
repair 


BCM - Lack of equipment, tools, or facilities 


BCM - Lack of technical skills necessary to complete 
repair 


BCM - Lack of parts 


BCM - Fails check and test, and maintenance is allowed 
to conduct check and test only 


BCM - Lack of technical data 

BCM - Bevond the authorized repair depth 
BCM - Administrative 

BCM --Condemned, repair not feasible 


[tems of repairable Material or Weapon: Support 
Svstem Discrepency Checked Mo Repair Required. This 
code 1s used for all discrepencies which are checked 

and found that either the reported deficiencv cannot 

be duplicated, or is operating within allowable 
tolerances. 


Repair or replacement of attaching units, seals, 
gaskets, etc., that are not integral parts of work 
unit coded itemis. 


This code is used when.a repairable item of material 
identified by a work unit code 1s repaired. 


This code 1s used to closeout a VIDS’MAF when .. 
component repair 1s to be performed at another facility. 


Failure of components/items undergoing test. 


This code is used when an item is calibrated and found 
serviceable without need for adjustment. 


This code is used when an item must be adjusted to 
meet calibration standards. 


This code is used when a maintenance action must 
be stopped or delaved while awaiting parts which are 
not available locally, and a component goes into 

an awalting parts status. 


This code is used bv an organizational activitv 

when it becomes nécessarv to closeout a maintenance 
action during or at the end of a reporting period 

Or anv reason. 


This code is entered when an item of material is 
removed and onlv the removal is to be accounted for. 
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This code is entered when an item is installed and 
only the installation action is to be accounted for. 


This code is entered when an item of material is 
removed due to suspected malfunction and the same or a 
like item is reinstalled. 


This code is entered when an item of material is 
removed to facilitate other maintenance and the same 
item is reinstalled. 


This code is used when an item of materialis . 
removed and replaced for a cannibalization action. 


This code is used when the time expended in locating 
a discrepency is great enough to warrant separating 
the troubleshoot time from the repair time. 


This code is used when actuallv treating corroded 
items, and includes cleaning, treatment, priming and 
painting. 


APPENDIX C 
WHEN DISCOVERED CODES 


DESCRIPTION 


Before flight - mission aborted - failure discovered bv aircrew 
Before flight - no abort - failure discovered by aircrew 
Inflight - mission aborted 

Inflight - no abort 

After flight/' between flight - failure discovered by aircrew 
Pilot/NFO weekly inspection 

Acceptance; transfer inspection 

Between flights - failure discovered by ground crew 
Daily inspection 

Preflight. daily, preflight, postflight, or turnaround inspection 
Special inspection 

Calender odd/major phase inspection 

Calender even inspection 

Administrative 

Functional checkflight 

Conditional inspection 

Quality assurance inspection 

Oil analvsis inspection 

Modification, SDLM, Overhaul' Airline mamtenance 
Related maintenance action 

In-shop repair; disassembly for maintenance 

Test bench/engine test stand operation 


Receipt or withdrawal from supplv 
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APPENDIX D 
MAXIMUM LIKELIHOOD DERIVATIONS 


l. MAXIMUM LIKELIHOOD ESTIMATORS 
a. Properties 
The method used in this thesis to estimate the model parameters is maximum 
likelihood. As such, it is necessary to review two important properties of this estimator. 
First, if a random sample of size n is taken of a random variable X whose probability 
mass function, p(x), depends on an unknown parameter, Û, then the probability law for 


the maximum likelihood estimator is asymptotically normal (as n gets large) with: 
MEAN = 0 and VARIANCE z (nK?y! (eqn D.1) 


where, K^» E[((0/00)Inp(x))"] for a discrete random variable X, and n equals the total | 
“number of observations. Thus, the maximum likelihood. estimator is asymptotically 
unbiased; that is, the expected value of the estimator is equal to the unknown 
parameter 0. Furthermore, the variance of the estimator is the Cramer-Rao lower 
bound, and represents the smallest possible variance of an unbiased estimator. 
INIGESTO: pp. 372,379) 

The second property describes the probability law for maximum likelihood 
estimators when the A $ for 12 l,..,n are not identically distributed. In the models 
discussed in this thesis, the probability mass function for each X, was dependent on 
One or more other random variables. Nevertheless, the distribution of the maximum 


likelihood estimator(s) 1s still asymptotically normal with: 


MEAN = @ and VARIANCE = (YK)! (eqn D.2) 
i 


where, gds E[((0/00)Inp(x))*] is evaluated using the specific parameters defining the 
ith observation. The difference in the variance in equation D.1 and the variance in 
equation D.2 is caused by probability mass functions for the X:s which are not 
identically distributed. For a complete derivation of this result, the reader is referred to 
W. J. Heintzelman's 1975 paper on Derermining the Failure Rate When Failure Times 


Are Not Known Exactly [Ref. 11]. 
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b. Confidence Intervals 
Confidence intervals for the unknown parameters can then be calculated using 
the maximum likelihood's probability law described in the previous section. These 


limits are defined as: 
CONFIDENCE LIMITS = 0 17 x 0g) 


where, 9’ is the maximum likelihood estimator for 9, Gg is the standard deviation of 
the maximum likelihood estimator, and 21-4/2 is the 1-4/2 quantile of the standard 
normal distribution. | 

c. Formulation of Probability Mass Functions 
The generalized format used to construct each model's probability mass function for X; 


1S: 
p(x.) = [P(failure occurred)"; x [P(right hand censored)]%i 


The random variable, X., is an indicator variable equal to zero if the data entry 
contains a failure, and equal to one if the observation is right band censored.? The 
following variables are used to derive the maximum likelihood estimate(s) in each of 


the models: 


° n : Total number of records in the final data set 
° m : Total number of records containing component failures 
e nm : Total number of records which were right hand censored 


2: EXPONENTIAL MODEL 
a. Properties 
The probability density function f(t), survival function S(t), and hazard rate 


function r(t) for the Exponential distribution are defined below. 
O) = Ae kt 1>0 
-()؟‎ ٢:) "لدم جح(‎ 


28D 


r(t) S(t) 


.. Richt hand censored observations identify those components that did not fail 
prior to the termination of the deplovment. 
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b. Formulation of the Maximum Likelihood Equation 
This section outlines the formulation of the maximum likelihood equation 
used to find A”, the maximum likelihood estimate for the failure rate parameter À. The 
final data set contains information supporting two record types. The first record tvpe 
identifies observations of components failing during some specified interval of time, T; 
to Tı +t. The probability equation representative of this type of record is given in 


equation D.3. 
P(T, St <T,+t,)= exp(-AT;) X (1-exp(-At,)) (eqn D.3) 


The second record type pertains to an observation of a component surviving 
past some time T;. These censored observations contain no information regarding the 
actual time of component failure. Equation D.4 defines the probability relationship 


associated with this type of observation. 
S(T.)= exp(-AT.) | (eqn D.4) 


The probability mass function, p(x;), for the Exponential model is obtained by 
combining equations D.3 and D.4. The random variable X, equals zero when a failure 
is observed and one when the data is right hand censored (i.e., component survived 


until the end of the deployment). 
p(x)= [exp(-AT)P x [exp(-AT;) x (1-exp(-Xt) 1-2 (eqn D.5) 


where, 
o ON) O; and T., 20 
e T; equals the aggregate flight time without component failure for the ith record 


e t; equals the flight duration of the sortie during which the failure was observed 


The hkehhood function for p(x.), L(p(x;)) can be expressed as the product of 
the n probability mass functions. The maximum likelihood estimate for A is the value 
of À” which maximizes L(p(x;)). Since the logarithmic function is a monotonically 
increasing function, the A’ which maximizes L(p(x;)) is identical to the A’ which 


maximizes the log-likelihood function K(p(x;)). The optimal value for A is found bv 
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taking the derivative of the log-likelihood function with respect to A, setting it equal to 
zero, and solving for A. 


Likelihood function: 


L(p(x,)) = F [exp(-AT,)I% x [exp(-AT;) x (1-exp(-At;)]* (eqn D.6) 
1 


17 1,2,3,...,n 
Log-likelihood function: 


K(p(x))* Y - ATx- AT(1-x) * (1-x)log[1-exp(-At,)] (eqn D.7) 
i ni 


The partial derivative of the log-likelihood function with respect to Å is: 


[( I-x.)texp(-A.t.)] 


00۳07 (eqn D.8) 


OK ` 
a IS Ix) 
1 
Equation D.9 is the result of setting Equation D.S equal to zero and 
rearranging the summations to eliminate the x;(s). The left hand side of Equation D.9 
is Summed over the set (1) of all records, censored and uncensored. But, the right hand 


side is only summed for the uncensored observations (j). 


texp(-At.) 77۲٤ 


» T, = > 20ھ‎ SS | (eqn D.9) 
: | l-exp(-At) 1=1,2 5 ۶۶۷۹ 


Cc. Maximum Likelihood Variance 
As stated in Section D.l, the maximum likehhood estimator for À is 
approximately normally distributed with mean equal to À and variance equal to bis 


for 17 1,...,n. [In this model: 


(t, exp(-At.) - S 
Isexp(- At») 


220 
RG 


1 


= x [exp(-XT;) X exp( -At;)] + T;i“exp(-À T.) 
where, - is evaluated using the maximum likelihood estimate of À (i.e, A=). The 
standard deviation of the maximum likelihood estimator and associated confidence 


limits displayed in Table 7 were computed using the equation listed above. 


3 WEIBULL MODEL 
a. Properties 
The probability density function f(t), survival function S(t), and hazard rate 
function r(t) for the Weibull distribution are displayed below. 


Rt) = 4 at®lexpf-(ay@) A, a>0 t20 
S(t) = exp[-(At)7] 
r(t) = Ade! 


b. Formulation-of the Maximum Likelihood Eguation 

This section describes the formulation of the eguation used to find the 
maximum likelihood estimates of the Weibull scale (A) and shape (4) parameters. The 
two different record types supporting this equation are identical to those identified in 
the preceding section. The key difference between the two models 1s that the Weibull 
X distribution possesses the “memoryless” property for a equal to one only. As such, the 
conditional survival probability is usually not equal to the unconditional survival 
probability. The conditional survival probabilitv for the Weibull distribution is defined 


as:10 


exp (AT + J 


"TC exp TA 


The conditional probability of failure during the interval T; to T. t, can be 
expressed as l-ST (t), the complement of the conditional survival probability. The 
probabilitv mass function for the Weibull model is then formed bv combining the 
conditional probability of a component failure within the interval T, to T, +t, (x, = 0) 


with the probability that a component survives time nb 
P(X.) = SA x expi(A T) Si (eqn D.11) 
1 


where, 
TON 0rand T, t 20 


e T; equals the aggregate flight time without component failure for the ith record 


!0The Subscripted T; is used to indicate the conditional survival probabilitv 
associated with a component that has survived T; hours of operation. 
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e t; equals the flight duration of the sortie during which the failure was observed 


The log-likelihood function, K(p(x;)), is: 


. 


K(p(x))= > (1-x)(log[1-S+ (t)]] + ((۸۲)-آ×‎ (eqn D.12) 
1 1 


+7777 


The partial derivatives of the log-likelihood function, K(p(x;)), with respect to the 
parameters À and d are displaved below. They have been simplified by removing the 
X;(s) and summing over the appropriate records. [he set 1= 1,2,3,...,n represents all data 


records (censored and uncensored), and the set J# m+ l,m F 2,...,n includes Ons 
records which are censored. 


> a-1 0 0 
OK ې‎ SA MT) Type) | Y اوی‎ 
oO" Lë i ٢ 
1 
٧ nd | ya | 
= _ y A E t)] P 70 log(AT;) Sy (t;) _ FAT loet 
00 |-ST(t) . J | | 
1 


c. Maximum Likelihood Variance 


The following equations were used to compute the variance for the maximum 
likelihood estimators, A and q”. These variances were then used to construct the 
confidence limits displayed in Table 8. 


2 : E 
K,“ with respect to À“: 


(IT, te) -T Sr 


KAA) = (ga 6h? x | 
i ()7( Jw ST) 


- 2 
TEN 


n er: 
K," with respect to 4: 


AICTE og T; * t) EO Ty" og. Ty) Sr (t)‏ د 


KG 
al I-ST (tj) | 


where, D equals KAT: og(X T))]^exp[-( ( 
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4.  GEOMETRIC MODEL 
a. Properties 
The probability mass function, Pr(F=s), and survival function, Pr(Fz s), for the 


Geometric distribution are given below. 


Pr(F=s)=pq9l  0<p<1,q=1-p and s=1,2,3,... 
Pr(F > s)= q3! 


where, F is the random variable representing the observation of the first failure, and q 
is the probability of the component surviving the sortie. Equations D.13 and D.14 
illustrate the Geometric distribution’s “memoryless” and constant failure rate properties 


respectively. 


Pr(F> a+b/F> a) = (q3 bygd = q? = Pr(F> b) (eqn D.13) 


Pr(F=s)/Pr(F z s) 7 (pa*lyq*! 2 p (eqn D.14) 


b. Formulation of the Maximum Likelihood Equation 
This section formulates the equation used to find the maximum likelihood 
estimate of the constant failure probability (p) associated with any given sortie. For 
the Geometric model, the two event types which can occur include components failing 


during the sth th 


sortie, or components failing during some unknown sortie after tHe s 
sortie. The probabilitv expression associated with the first tvpe of event is given in 
equation D.15, and equation D.16 defines the relationship for those observations where 


the component survived at least s sorties. 


Pr(F=s) = وم‎ (EM DIS) 


Pr(F>s) = q? (eqn D.16) 


The probability mass function for this model is obtained by combining 


equations D.15 and D.16. Again, X is allowed to take the values zero and one only, 


with zero indicating a component failure on the sth 


th 


sortie and a one indicating a 


component failure after (be si" sortie, 
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p(x;) = [p(1-p)S 1] ^i x ((1-p)5]S (eqn D.17) 


where, 
° OSpSl, E b m 
e s: equals the aggregate number of sorties flown inclusive of the sortie causing 
the failure when x:=0 and equals the aggregate number of sorties without à 
component failure when x;= 1 
The maximum likelihood estimate p' of the unknown sortie failure rate will be 
found by taking the derivative of the log-likelihood function [K(p(x;))] with respect to 
p, setting the derivative equal to zero and solving for p. 


K(p(x))= 2 (1-x;)lIn(p) + (s;-1)(1-x;)In(I-p) t s;xin(l-p) (eqn D.18) 


OR (lx) [(s;-1)(1-x;)] Käre 
i P (1-p) (1-p) 


The partial derivative listed above can be simplified by removing the indicator 
variables and adjusting the summations to produce the closed form expression for p' 
given in equation D.19. 


Y =p El e 
J 1 


1 1 
ls l-p 





p = [T E | (eqn D.19) 


-14 


where, j= 1,2,3,....m 1s the set of uncensored records (1.e., component failures ES 
1= 1,2,3,....n 1s the set of all records, censored and uncensored. Thus, the maximum 
likelihood estimator p’ equals the number of flights experiencing a component failure 
divided by the total number of sorties flown. 
c. Maximum Likelihood Variance 

The maximum likelihood estimator’s variance was computed using the 
formulas described in Section D.l. The equation defining the Geometric model's 
expression for K;^ is given below: 
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KZ = (1-1 p + spy? (egn D.20) 


where, ۴2 is evaluated at p=p’. These variances were used to compute the confidence 


limits listed in Table 9. 


5. MIXED MODEL 

a. Properties 
Combining a constant sortie failure probability, p,, with a flight hour probability 
distribution results in an expression for the probability of failure, p,, for the ¡th sortie. 
This is expressed as the complement of a component surviving both the stresses of the 


sortie and the flight hours flown: 
p; — l-expl(-At.)(l-p,) —9'$po $1, &»0,t 20 (eqn D.21) 


where, the flight hour failure distribution is Exponential. 

The model is then formulated in the same manner as the Geometric model 
described in the preceding section. The probabilitv that the first failure occurs on the 
sth sortie is expressed in equation D.22. Note that the Pi(s) are not necessarilv equal 


for 1= 1....,s-1. 


(eqn D.22)‏ ا 2 I (epy)‏ 7 وس ات 


Equation D.23 results from using the expression for p; given in equation D.21 


to rewrite equation D.22. 
Pr(F=s) = [1L-exp(ADA-P Map) Pexp(-AT)] dl EDS) 


b. Formulation of the Maximum Likelihood Equation 
This section describes the formulation of the equation used to find the 
maximum likelihood estimates for the failure rate parameter (A) and the constant sortie 
probability of failure (pg). The probability relationship supporting the observation of 
the first failure on the sth sortie is given in equation D.23, and the expression for the 


th 


observation of a failure at some unknown time after s" sortie is listed below. 
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Pr(F>s) = II (1-p¡) = exp(-AT)1-po)° (eqn D.24) 
i 


The probability mass function p(x;) for this model results from combining 
equations D.23 and D.24. Notice that when À equals zero the probability mass 
function for the Mixed model is identical to the probability mass function expressed for 
the Geometric model in equation D.17, and when p, equals zero the probability mass 
function is equivalent to the. Exponential probability mass function given in equation 
D.5. The above statement implies that the Mixed model has the capability of 


modelling the Geometric or Exponential cases as the situation would require. 


e 


p(x)= [expte Ma? Y Tex (AT INS x [(1-p, Si exp(-AT FS 


where, 
° 051 ٢ (٢ | s.= 1,2,3,… and T., t 20 
1 1I 
e T; equals the aggregate flight time without component failure for the ith record 
e t equals the flight duration of the sortie during which the failure was observed 


° - equals the aggregate number of sorties flown inclusive of the sortie causing 


5 
the failure when x:=0 and equals the aggregate number of sorties without à 
component failure When x; = 1. 

The log-likelihood function K(p(x;)), and the partial derivatives of the log- 
likelihood functions with respect to the parameters À and Po are given below. The 
maximum likelihood estimates for À and p, were obtained utilizing the Quasi- Newton 


method described in Appendix G. The log-likelihood function, K(p(x;)). is: 


K(p(x;))= >( 1-x;)(In[1-exp(-At.)X1-p.)] +[(s-1)in(I-p.)]-K Tj +x;([s.In(1-p.)]-A T; 
1 


17 1.2,3,....n 


The partial derivatives of the log-likelihood function, K(p(x.)), have been simplified by 
removing the x,(s) and summing over the appropriate records. [he set 1— 1,2.3.....n 
represents all data records (censored and uncensored), and the set j=m+ a 


includes only those records which are censored. 


10 





OK 1-p_)t;exp(-At, 
ŠK, (poletel) جح‎ (eqn A.25) 
On ¡ l-exp(-At; )( ISP, j J 
0 -At; .-1 
SS 9۹ (eqn A.26) 


Opo 7 lexp-AtXl-po) l-Po j Po 


c. Maximum Likelihood Variance 
The equation for the variance of the probability distribution for Po IS 
identical to equation D.20 when À equals zero. The equation for the K. for À is given 


below. 
> 2 r i e . 
A =I 01 TA x po(1-po) i ېا‎ E )5۔1‎ 


where, Kea is evaluated at A=A'=0. These variances were used to compute the 


confidence limits for. À listed in Table 9. 


p 


sal 


APPENDIX E 
SINGLE PARAMETER OPTIMIZER 


COMPUTES MAXIMUM LIKELIHOOD VALUE FOR A SINGLE UNKNOWN 
PARAMETER USING GOLDEN RATIO LINE SEARCH. 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/RL/FTNF , FTF, PARAM, PSD 

COMMON/ INT/ CENSOR , NF, COUNT 

DIMENSION CENSOR (2009 ,NF (200) ,FTNF (200) ,FTF (200) ,PARAM (2) 
INTEGER CENSOR, COUNT, TFLAG 

REALX8 LB,LC 


DATA LB/.0001D0/ITER/1/K/1/STEP/.00001D0/ALPHA/.0001D0/ 


INITIALIZE 
COUNT=1 
PARAM(1)=. 00 
TFLAG=0 


COMPUTATION OF SEARCH INTERVAL TAU. TAU IS APPROXIMATELY 
.618 AND TAU] IS APPROXIMATELY .382 

TAUZ(5**,5-1)/2 

TAU1=1-TAU : 


INPUT DATA 
READ (3,50,END=99) CENSOR (COUNT) ,NF (COUNT) , FINF (COUNT) , FTF (COUNT) 
COUNT=COUNT+1 
GO TO 5 
COUNT=COUNT-1 


BOUNDING PHASE - SWANN'S METHOD 
CALL ML(LB F 
ALPHA= ALEHA (2 رات‎ 
CALL ML(ALPHA,F2 


IF NEW VALUE IS LESS THAN PREVIOUS VALUE TERMINATE BOUNDING 
EE (7 1 ٥ 
F1=F2 
K=K+1 
GO TO 6 
UB=ALPHA 
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*DEFINE INITIAL POINTS & ASSOCIATED FUNCTIONAL VALUES 


B1=TAU1*(UB-LB)+LB 
B2=TAU*(UB-LB)+LB 
6۸8۱۲۱٢٢۱٤4 
CALL ML B2,F2 


COMPARES THE FUNCTIONAL VALUES OF THE INTERNAL PTS. Bl & B2, AND 
SELECTS THE VALUE WITH THE LARGEST FUNCTIONAL VALUE AS THE UPDATED 
PARAMETER ESTIMATE 


IF (POT: Li THEN 
PARAM(2 
UB=B2 
82-1 
F2=F1 
B1=TAU1*(UB-LB)+LB 
CALL ML(B1,F1) 
ITER=ITER+1 
GO TO 200 

END IF 


2 


A A 


150 


MES) 


TE an GE.F1) THEN 
PARAM(2)=B2 
LB=B1 
Bl-B2 
F1=F2 
o p eg 
٤٢ ٩272 
٣۲85۲1 
GO TO 200 
END IF 


200 CALL TERM(F1,F2,TFLAG) 
PARAM(1)=PARAM(2) 
796۶۱6٣۲۸6 ۰0۰7۰۰ TO 275 
GO TO 150 


* WRITE OUTPUT 
275 WRITE (6,300) 
WRITE 9 310) PARAM(1) 
CALL SD 
WRITE (6,305) PSD 
leč paie. .96*PSD 


LC=PARAM(1)-1.96*PSD 
WRITE (6,320) LC,UC 
^ FORMAT 


50 BORSE (IM IX I4, 5X,F6.1,1X,F3.1) 
300 FORMAT ('OWUC: ' / OCONVERGENCE CRITERION ESTABLISHED!) 
305 FORMAT ('OSTANDARD DEVIATION IS | = EEG 
310 FORMAT ('OMLE FOR UNKNOWN PARAMETER IS: ' F9.8 
320 FORMAT" 0.95 ASSYMPTOTIC CONFIDENCE LIMITS ARE :',F5.4, 
2999 S l 

N 


SUBROUTINE TERM(F1,F2,TFLAG) 

* DETERMINES IF TERMINATION CONDITION HAS BEEN SATISFIED 
IMPLICIT REAL*8 (A-H,0-Z) 
COMMON/RL/FTNF , FTF, PARAM, PSD 
COMMON/INT/ CENSOR , NF, COUNT 
DIMENSION CENSOR(200) ,NF(200) , FTNF(200) , FTF(200) , PARAM(2) 
INTEGER TFLAG 


TEST=ABS (F1-F2) 
IF (TEST.LE..00003) TFLAG=1 


RETURN 
END 
SUBROUTINE ML(PT,FVAL) 

* COMPUTES LOG-LIKELIHOOD FUNCTIONAL VALUE 
IMPLICIT REAL^8 (A-H,O-Z) 
COMMON/RL/FTNF ,FIF,PARAM E 


COMMON/INT/CENSOR, NF , COUN 
DIMENSION CENSOR(200) ¿NE (200), FTNF(200),FTF(200), PARAM(2) 


INTEGER CENSOR, COUNT, TFLA 
REAL%8 Kl 
FVAL=0 


* INPUT FLIGHT DATA 
DO 5 I=1 COUNT 


x UNCENSORED RECORDS 


ZB 


IF (CENSOR(I).EQ.0) THEN O. 
K1=EXP(-PT*FIF(I)) 
FVAL=FVAL+LOG(1-K1)-PT*FTNF(I) 
GO To 5 

END IF 


CENSORED RECORDS 
FVAL=FVAL-PTXFTNF(I) 
CONTINUE 


RETURN 
END 


SUBROUTINE SD 


COLLO ASSYMPTOTIC STANDARD DEVIATION FOR THE MAXIMUM LIKELIHOOD 
511 


IMPLICIT REAL*8 (A-H,O-Z) 

COMMON/RL/FTNF,FTF,PARAM,PSD 

COMMON/INT/CENSOR,NF,COUNT 

DIMENSION CENSOR(200) ,NF(200),FTNF(200),FTF(200) , PARAM(2) 
INTEGER CENSOR, COUNT 

8 +4 +1 9 


SUMKI=0. 


INPUT FLIGHI DATA 
DO 5 I=1,COUNT 


UNCENSORED RECORDS 

IF (CENSOR(I).EQ.0) THEN 
M RETE CT) 
K2=EXP(-PARAM(1)*FINF(1)) 
A TEA 
KISO=K3*(K2*(1-K1))+(FTNF(I)* 
Go TO 6 

END IF 


CENSORED RECORDS 
KISO=FTINF(I)*FINF(I)*K2 
SUMKI=SUMKI+KISO 
CONTINUE 


PSD=SQRT(1/SUMKI) 


RETURN 
END 


x 


)) ({ETE(I)*K1/(1-K1))-FINF(1) 
RENE 


SUBROUTINE SDI 


COMPUTES ASSYMPTOTIC STANDARD DEVIATION FOR THE MAXIMUM LIKELIHOOD 
ESTIMATE ہ0‎ ]8٦. 505۰ 6 ۶ +077۶7 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/RL/FTNF , FIF , PARAM, PSD 

COMMON/ INT/ CENSOR ,NF , COUNT 

DIMENSION CENSOR(200$ ,NF(200) , FTNF(200) , FTF(200) , PARAM(2) 
INTEGER CENSOR, COUNT 

REAL*8 L,LL1,LL2 


INITIALIZE PARAMETER ESTIMATES 
P=.0232D0 
L=.0000D0 
SUMKI1=0.0DO 
SUMKI2=0.0DO 
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x INPUT FLIGHT DATA 
DO 5 I-1,COUNT 
Ci= ga P*NE(I)) 


C2=(1- P)*^NFI) 
63 - 71 
C4=C3/(1- P 


C5=C4/(1-P 

Cez (^ P)*FTF(I))/P- AE 

501101۰01111+ )617 0(0 ا د مرگ‎ LA 

SUMKI2=SUMKI2+C6*C6*P*C3+FTNF(1)*FTNF(1)*C2 
CONTINUE 


COMPUTE STANDARD DEVIATIONS 
PSD1= DE 
PSD2-SORT(1/SUMKI2 


7 COMPUTE CONFIDENCE LIMITS 
LL1=P-1.96*PSD1 
UL1=P+1.96*PSD1 
LL2=L-1.96*PSD2 
UL2=L+1.96*PSD2 


PUT 

WRITE Z LECL, UCI 
320 i DOS ۹5511۳۱) ٢ ٩ DENCE LIMITS FOR P: ',FS.4,',',F5.4 
C ۰ 

WRITE 1 ۹-7 CCZ UEZ 
330 Sd ' 0.95 ASSYMPTOTIC CONFIDENCE LIMITS FOR LAMBDA: ',F5.4,! 


RETURN 
END 


o 


HRS 


TIME BETWEEN FAILURES 


1000 


100 
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Exponential Quantile-Quantile Plot: WUC 46X 1600. 


APPENDIX F 
QUANTILE-QUANTILE PLOTS 
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Figure F.12 Weibull Quantile-Quantile Plot: WUC 7445M00, 
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Weibull Quantile-Quantile Plot: WUC 6918100. 


Figure F.14 
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APPENDIX G 
DUAL PARAMETER OPTIMIZER 


QUASI-NEWTON METHOD‏ .ا 

The Quasi-Newton method was used to determine the maximum likelihood 
estimates for the Weibull and Mixed models. This method was derived from the 
Newton's method through the Taylor's expansion of the objective function, f(x), where 
x is the nX 1 column vector containing values for the n parameters being estimated. 
The Newton method uses the quadratic approximation to Kx), equation G.l, evaluated 


at x) to iteratively determine an optimal solution. 
x) ~ Rx + Taxi) Tax + ہو لظکم 0۱2۸ا ود‎ (eqn G.1) 


The next point in the sequence is found bv forcing the gradient of the quadratic 


approximation (equation G.2) to zero and solving for x(K+ 1) 


Vi(x)= Vx + Vx) Ax = 0 | (egn G.2) 


Allowing Ax to egual x(K+ l). (K) equation G.2 can be solved for x(k + 1) 


x(k + D= y(K)_g 297K) yy (by (eqn G.3) 


The two-parameter optimizer displaved in Appendix G did not use the Newton 
method discussed above because the technique can be unreliable for nonquadratic 
functions. Instead, a more robust Quasi-Newton method was used. The major 
difference between the two methods is that the Quasi-Newton does not require use of 
the inverse Hessian matrix vin). Instead, it uses an nXn “A” matrix to estimate 


the inverse Hessian matrix. The AS) 


matrix was iteratively updated using the Brovden- 
Fletcher-Shanno method, and should be a fair approximation to ræv! when x* 1s 
the optimal solution. Substituting ACK) for ۹ (١ in equation G.3 results in the 
Quasi-Newton method for the determination of the next point in the sequence. This 
iterative process continues until the objective function is maximized. [Ref. 12: pp. 


112-114] 


90 


FORTRAN PROGRAM 


COMPUTES MAXIMUM LIKELIHOOD ESTIMATES FOR TWO PARAMETER 
MODELS USING QUASI-NEWTON METHOD. 


IMPLICIT REAL*8 (A-H,0-Z) 
COMMON/RL/ALPHA,FVAL,UFVAL,DIR,FTNF,FTF,GRAD,PARAM 
COMMON/INT/CONFLG, COUNT, TFLAG, CENSOR,NF 

DIMENSION CENSOR(200),NF(200),DIR(2), 'ETNF(200) 
DIMENSION FTF(200), ,GRÀD (4) ,PARAM(4 

DIMENSION A(4),A1(4),A2 DELP(2),DELG(2) 

INTEGER CENSOR, CONFLG, COUNT, TFLÀG 


INITIALIZE 'A' - APPROXIMATION FOR INVERSE HESSI 
DATA A(1)/1.0D0/A(2)/0.0D0/A(3)/0. 050/A(2)/1. ODO/ITER/0/IRES/21/ 


CONFLG=0 

DIR 0." 

DIR(2)=-1.0D0 
INPUT DATA 

CALL INPUT 


WRITE (6,*) 'COUNT' COUNT 


OUTPUT STARTING VALUES 
WRITE 0 
WRITE (6,500) (PARAM(I),I=1,2) 


COMPUTE LOGLIKELIHOOD VALUE FOR INITIAL STARTING ESTIMATES 
CACC LLIKE(1) 
FV=-1*FVAL 
WRME( 6,510) FV 


CHECK FOR CONVERGENCE 
CALL CONCHK(1) 


IF (CONFLG.EQ.1) GO TO 1000 


ITER=ITER+1 
WRITE (6,520) ITER 


IF ITERATIONS EQUAL FIFTY, PROGRAM IS TERMINATED EARLY 
TF (ITE GE.50) 1T 
WRITE (6, S60). 
8 9 ۵۷ 
END IF 


IF IRES=1, OBJECTIVE FUNCTION IS NO LONGER IMPROVING AND GRADIENT 
RR OL COULD IO ZERO. THE "A" MATRIX IS RESET TO THE IDENTITY MATRIX. 
PILE BESKØD > TRES) THEN 
A(1)=1.0D 


A(2)=0.0DO 

A(3)=0.0D0 

A(4)=1.0DO 

IRES=IRES+20 
END IF 


BOSUSPOTE DIRECTION OF IMPROVEMENT 
md DECR EM DECRE AM 


DIR(2)=-1*(A(3)*GRAD(1)+A(4)*GRAD 
COMPUTE LINE STEP SIZE (ALPHA) 
CALL LINE 


COMPUTE DIFFERENCE BETWEEN PARAM(K+1) AND PARAM(K), AND 
DIFFERENCE BETWEEN GRAD(K+1) AND GRAD(K) 
000 1, 2 
Ne: =PARAM(I+2)-PARAM(I) 
DELG(I)=GRAD(I+2)-GRAD(I) 


91 


** N 


CONTINUE 


10 


* UPDATE 'A' MATRIX 


C1=DELP(1)*DELG(1)+DELP(2)*DELG(2) 


A1(1)=DELP(1)*DELP(1)/C1 
A1(2)=DELP(1)*DELP(2)/C1 
A1(3)=DELP(1)*DELP(2)/C1 
A1 (4)=DELP(2)*DELE cl 


C2=1- DEL Ds C1) 
C3=-1 EE 7ه‎ En 
C4=-1*((DELP(2)*DELG C1 
C5=1 = ( (DELP (2 J DELG ٥٦ 
A2(1)=C2*(A(1)*C2+A(2)*XC3)+C3*(A(3)*C2+A(4)*C3 
A2(2)=C2*(A(1)*C4+A(2)*C5)+C3*(A(3)*C4+A(4)*C5 
A2(3)=C4*(A(1)*C2+A(2)*C3)+C5*(A(3)*C2+A(4)*C3 
20ھ‎ 6 ٦ AC4+A 2)XC5)+C5*(A(3)*C4+A(4)*C5 
DO 20 I=1 
A(I) 2x2(1)+A1(1) 
CONTINUE 


20 


* UPDATE OLD PARAMETER ESTIMATES AND GRADIENTS TO NEW PARAMETER 


* ESTIMATES AND GRADIENTS 


DO 25 I=1,2 
PARAM(I)<PARAM(It2) 
GRAD(I)=GRAD(I+2) 

CONTINUE 


25 


* OUTPUT UPDATED PARAMETERS 


WRITE (6,500) (PARAM(I),1=1,2) 


* OUTPUT FUNCTIONAL, VALUE 


FV=-1*FVAL 
WRITE(6,510) FV 


* CHECK FOR CONVERGENCE 


CALL CONCHK(1) 
IF (CONFLG.EQ.1) GO TO 1000 


* START NEW ITERATION 
GO TOFI 


* PROGRAM TERMINATION 
1000 WRITE (6,530) 


WRITE 6, E 
WRITE (6,550) ñ 


* OUTPUT FORMAT 


FORMAT ('OSTARTING VALUES:',/,' ---------------- 
FORMAT ('OPHAT IS: ',F15.8,/,' LHAT IS: ',F15.8) 
FORMAT ('OLOGLIKELIHOOD VALUE IS: pu 
FORMAT ('OITERATION: ',I13,/,0 2222222355 | 

FORMAT ('OPROGRAM COMPLETION') 

FORMAT ('OFINAL "A" MATRIX IS:' 

FORMAT ('0',2(F15.8),/,' ' 


OEIL SD 
FORMAT ' OUNSÜCCESSFÜL TERMINATION. ITERATIONS = 50') 


STOP 
END 


SUBROUTINE INPUT 


490 
500 
510 
520 
530 
540 
550 
560 


READS FLIGHT AND MAINTENANCE RECORDS USED IN THE CALCULATION 
OF MIXED MODEL MAXIMUM LIKELIHOOD ESTIMATES FOR PHAT AND 
LAMBDA HAT 


IMPLICIT REAL*8 (A-H,0-Z) 
COMMON/RL/ALPHA,FVAL,UFVAL,DIR, Ge FTF,GRAD, PARAM 
COMMON/INT/CONFLG, COUNT, TFLAG, CENSO R,NE 

DIMENSION CENSOR(200) ,NF(200), IDIR(2), ,FTNF (200) 
DIMENSION FTF(200), GRAD(4), PARAM(4 ) 

[INTEGER CENSOR, CONFLG, COUNT, TFLAG 


* STARTING POINTS FOR PARAM(1),PHAT, PARAM(2),LHAT, AND RECORD COUNTER 
PARAM (1 = Ro 


COUNT=1 


* INPUT FLIGHT DATA 
5 na (3,500,END=1000) CENSOR(COUNT),NF(COUNT),FTNF(COUNT),FTF(COUN 
T 


COUNT=COUNT+1 
GO TO 5 


x INPUT FORMAT 
500 HESRMAT O (11,IX/14,5X,F6.1,1X,F3-1) 


1000 COUNT=COUNT-1 
RETURN 
END 


x xx 


SUBROUTINE CONCHK(N) 
* DETERMINES IF OPTIMAL SOLUTION HAS BEEN FOUND 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/RL/ALPHA, FVAL,UFVAL,DIR,FTNF,FTF,GRAD, PARAM 
COMMON/INT/CONFLG, COUNT, TFLAG, CENSOR,NF 

DIMENSION CENSOR(200),NF(200), ,DIR(2), FTNF (200) 
DIMENSION FTF(200),GRAD(4), PARAM(4) 

INTEGER CENSOR, CONFLG, COUNT, TFLAG 

REAL*8 NORMG 


k TERMINATE IF THE NORM OF THE GRADIENT (NORMG) IS LESS THAN THE 
* THE SPECIFIED CRITICAL VALUE (CRIT) 
NORMG- ( (GRAD (N) *GRAD(N) )- (GRAD(N*1)*GRAD(N*1)))**.5 
IF (NORMG.LE.1.) THEN 
CONFLG=1 
WRITE (6,500). 
D IF 


FTEST=ABS (UFVAL-FVAL) 


IF (FTEST.LE..000001.AND.NORMG.GT.1.) THEN 
CONFLG=1 
WRITE (6,510) 

END IF 


* OUTPUT FORMAT 
500 FORMAT TIONS RMAT TERMINATION - OPTIMAL FOUND!) 
S10 FORMAT ('OPROGRAM TERMINATION ~ LACK OF FUNCTIONAL IMPROVEMENT ' ) 


RETURN 
END 


SUBROUTINE LINE 
7 COMPUTES MAXIMUM STEP SIZE (ALPHA) 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/RL/ALPHA, FVAL,UFVAL,DIR,FTNF,FTF,GRAD, PARAM 
COMMON/INT/CONFLG, COUNT, TFLAG, CENSOR,NF 

DIMENSION CENSOR(200),NF(200), ,DIR(2); ,FTNF (200) 
DIMENSION FTF(200) ,GRAD(4) , PARAM(4) 

INTEGER CENSOR, CONFLG, COUNT, TFLAG 

REAL*8 LB 


UFVAL=FVAL 
K=1 


TFLAG=0 
ALPHA=0 .0DO 
ALPHAK-.0001D0 
LB=0.0DO 


x COMPUTATION OF SEARCH INTERVAL TAU. TAU IS APPROXIMATELY 
* .618 AND TAU1 IS APPROXIMATELY .382 

TAUZ(5**.5-1)/2 

TAU1=1-'TAU 


č BOUNDING PHASE - SWANN'S METHOD 
STEP=.00001DO 
RE لل‎ 


PARAM(4)=PARAM(2)-STEP*DIR 
CALL LLIRE(3) 
FI=FVAL 
PARAM DOMI E 
PARAM(4)=PARAM(2)+STEP*DIR(2 
CALL ILIKE) 
F2=FVAL 

IF (Fl.GE.F2) THEN 
F1=F2 
GO TO 5 

END IF 

IF (F2.GT.F1) STEP=-1*STEP 

5 ALPHA=ALPHA+ (2**K) *STEP 
PA s. a 
PARAM(4)=PARAM(2.)+ALPHA*DIR(2 


k — COMPUTE LOGLIKELIHOOD VALUE 
CALL LEIKE (8) 
F2=FVAL 


< IF NEW VALUE IS LESS THAN PREVIOUS VALUE TERHINATE BOUNDING 


IF (F2.GT.Fl1) GO TO 10 
Fl<F2 

K=K+1 

GO TO 5 


10 UB=ALPHA 
x CONDUCT GOLDEN SECTION SEARCH 


B1=TAU1*(UB-LB)+LB 

E EL‏ خر 
PARAM(4)=PARAM(2)+B1*DIR(2‏ 
CALL LLIKE(3)‏ 

FB1=FVAL 

B2=TAU*(UB-LB)+LB 

PARAM کہ‎ +5 
PARAM(4)=PARAM(2)+B2*DIR(2 
CALL LLIKE (3) 
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FB2=FVAL 
50 IF (FB1.LT.FB2) THEN 

ALPHA=B1 
UB=B2 
82-81 
۶82-81 
B1=TAU1*(UB-LB)+LB 
Baal BE ٦ 
PARAM(4)=PARAM(2)+B1*DIR(2 
CALL LLIKE(3) 
FB1=FVAL 
GO TO 100 


END IF 
IF (FB2.LE.FB1) THEN 
ALPHA=B2 
LB=B1 
B1=B2 
FB1=FB2 
B2=TAU*(UB-LB)+LB 
E RAE دص‎ 
PARAM(4)=PARAM(2)+B2*DIR(2 
PARE LLIKE(3) 
FB2=FVAL 
GO TO 100 
END IF 


* CHECK FOR LINE SEARCH TERMINATION 
100 CALL TERM(ALPHAK) 
ALPHAK=ALPHA 
Pole LAG. 50.1) GO TO 5000 
CO TOSO = 


5000 PARAM(3)- PARAM(1)+ALPHA*DIR(1) 
PARAM(4)=PARAM(2)+ALPHA*DIR(2) 

CALE LLIKE(3) 

RETURN 

END 

SUBROUTINE TERM(ALPHAK) 
* DETERMINES IF CONVERGENCE CRITERIA IS ESTABLISHED FOR LINE 
* SEARCH TERMINATION 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/RL/ALPHA, FVAL,UFVAL,DIR,FTNF,FTF,GRAD, PARAM 
COMMON/INT/CONFLG, COUNT, TFLAG, CENSOR,NF 
DIMENSION CENSOR(200),NF(200), »DIR(2), 'FTNF(200) 
DIMENSION FTF(200),GRAD(4),PARAM(4) 

INTEGER CENSOR, CONFLG, COUNT, TFLAG 

TFLAG=0 

CTEST=ABS((ALPHAK-ALPHA)/ALPHA) 

IF (CTEST.LE..003) TFLAG=1 

RETURN 

END 

SUBROUTINE LLIKE(N) 


* COMPUTES LOGLIKELIHOOD FUNCTIONAL VALUE AND GRADIENTS FOR 
* FOR THE WEIBULL MODEL 


IMPLICIT REAL^8 (A-H,0-Z) 


COMMON/RL/ALPHA,FVAL,UFVAL,DIR,FTNF,FTF,GRAD,PARAM 
COMMON/ INT/CONFLG, COUNT, TFLAG , CENSOR , NF 

DIMENSION CENSOR(200) ,NF(200), »DIR(2), ,FTNF (200) 
DIMENSION FTF (200) ,GRAD(4) ,PARAM(4) 

INTEGER CENSOR, CONFLG, COUNT, TFLAG 


* ESTABLISH BOUNDARY CONDITIONS 
* PARAM(N)=SHAPE PARAMETER;PARAM(N+1)=SCALE PARAMTER 


TE Ea «LT. .00000000 1) PARAM(N)=.0001 
IF ۰)۲ Dr .0000000001) PARAM(N+1)=.0001 


FVAL-O. 
GRAD(N)=0 
GRAD(N+1)= st, 


% INPUT FLIGHT DATA 


DO 5 1 0۸7 


C1=PARAM(N+1)*FTNF(I) 
C2=C1**PARAM(N) 
٠۰) 9 70 


* UNCENSORED RECORDS 


IF 0" m THEN 
C3=PARAM(N+1)*(FINF(I)+FTF(I)) 
C4=C3**PARAM(N) 
CS=EXP(-C2 
C6=EXP(-C4 
FUNCTIONAL VALUE CALCULATIONS 
FVAL=FVAL-LOG(C5-C6) ° 
GRADIENT CALCULATIONS 
7ء‎ Ee EE BE (CEXC4*LOG(C3)))/(CS 
GRAD(N+1) LESEN **PARAM(N))-C6 E 
E ) PSPARAM (N) ) )*( PARAM(N)*PARAM(N*1)**(PARAM(N)-1)))/(C 
TO 


END TF 


759) +FTF (I 
E 


CENSORED RECORDS 


FUNCTIONAL VALUE CALCULATIONS 
FVAL=FVAL+C2 

GRADIENT CALCULATIONS 
GERD UDE =GRAD(N)+(C2*LOG ( Eu 
GRAD(N+1)=PARAM(N)*(FTNF(I)**PARAM(N))*(PARAM(N+1)**(PARAM(N)-1 


8 

CONTINUE 
RETURN 
END 


SUBROUTINE SD 


COMPUTES THE STANDARD DEVIATION FOR THE MAXIMUM LIKELIHOOD 
ESTIMATORS FOR THE WEIBULL DISTRIBUTION, WHICH JS USED 
TO CONSTRUCT .95 CONFIDENCE INTERVALS FOR THE SHAPE AND 
SCALE PARAMETERS. 


IMPLICIT REAL*8 (A-H,O-Z) 
COMMON/RL/ALPHA,FVAL,UFVAL,DIR,FTNF,FTF,GRAD,PARAM 
COMMON/INT/CONFLG, COUNT, TFLAG, CENSOR, NF 

DIMENSION CENSOR(200),NF(200), DIR(2), FINF(200) 
DIMENSION FTF (200) ,GRAD (4) , PARAM (4) 

INTEGER CENSOR, CONFLG, COUNT, TFLAG 

REAL*8 LL1,LL2 


SUMKI1=0. 


SUMKI2=0. 
C=(PARAM(1)*(PARAM(2)**(PARAM(1)-1)))**2 
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x + 


* * A *F 


x INPUT FLIGHT DATA 
DO 5 I=1,COUNT 


IF (FTF(I).EQ.0.) FTF(I)=.000001 
IF (FTNF(I).ÉQ.O. 1 FINF(I)=.000001 
CI=(FTNF(I FETEI ) 

C2=PARAM(2)*C1 
C3=PARAM(2)*FTNF(I) 
NEM í ) 
SUMKI1=SUMKI1+CX((((FTNF( 


P(-1*(C3**PARAM(1))) 
C HE 
5 


*PARAM(1)-C1**P EE TEN‏ و 
2س 
M(1‏ 
) 


7 -C4)+(FTNF(I)**PARAM(1))*(FINF(I 
C — )**PARAM(1))*EXP(-1^(C 1 


C 
kP 
I 


SUMKI2= تہ‎ CAM AM LOG c (C2**PARAM(1)*LOG(C2)))*C4) 
C. **2)/(1-C4)*(C3**PARAM(1)*LOG(C3)) 
5 CONTINUE 
|| اور ہے‎ 
PSD2=SQRT (1/SUMKI2 


A CONFIDENCE INTERVALS 
‘ LLI=PARAM(1)-1.96*PSD1 
UL1=PARAM(1)+1.96*PSD1 
FE 2=PARAM( 2)=1;96XPSD1 - 
UL2=PARAM(2)+1.96*PSD1 


A OUTPUT 
WRITE DIO ۲7 
RITE GE 110) PSD2,LL2,UL2 


100 ٣٥۷۰۰٦٦٦٠٠٠٠۶5 01:۰ ! F6.5,' CONEIDENCE INTERVAL ہ۶9"‎ ۶ ٠۶٦ 
10 A TEDL: 6٢و‎ ٢١ CONFIDENCE INTERVAL (',E6.5,',',F6.5,')' 
RETURN 
END 


SUBROUTINE LLIKE(N) 


COMPUTES LOGLIKELIHOOD FUNCTIONAL VALUE AND GRADIENTS FOR 
MIME MIXED MODEL. 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/RL/ALPHA, FVAL,UFVAL,DIR, FINF,FTF,GRAD, PARAM 
COMMON/INT/CONFLG, COUNT, TFLAG, CENSOR, NF 

DIMENSION CENSOR(200), NF (200), ik , FINF(200) 
DIMENSION FTF(200), GRAD(4), PARAM (4) 
1017000000505۰ CONFLG, COUNT, TFLAG 


* ESTABLISH BOUNDARY CONDITIONS 
IF (PARAM k nm PARAM(N)=.0001 
IF (PARAM(N).GE.1.) PARAM(N)=. 3 
IF (PARAM(N+1).LT..0000000001) PARAM(N+1)=.0001 


FVAL-O. 

GRAD (N) = =0. 
GRAD(N+1)=0 
Cl=1- PARAM(N) 
C2=LOG(C1) 


* INPUT FLIGHT DATA 
00757721 COUNT 


* UNCENSORED RECORDS 
IF (CENSOR(I). وت‎ THEN 


X f 








| C3=EXP (-PARAM(N+1)*FTF(I)) 
پر‎ FUNCTIONAL VALUE CALCULATIONS 
FVAL=FVAL-LOG(1-C3*C1)-(NF(1)-1)*C2+PARAM(N+1)*FTNF(I) 
x GRADIENT CALCULATIONS 


GRAD(N)=GRAD(N)-(C3/(1-C1*C3 F i 
کت 2ک , 2ء‎ 


O 


GO TO 5 
END IF 


* CENSORED RECORDS 
d FUNCTIONAL VALUE CALCULATIONS 
FVAL=FVAL-NF(I)*C2+PARAM(N+1)*FTNF(I) 
* GRADIENT CALCULATIONS 
موہ‎ s Ed 
GRAD(N+1)=GRAD(N+1)+FTNF(I) 
5 CONTINUE 


RETURN 
END 


98 


APPENDIX H 
WEIBULL SIMULATION OF COMPONENT FAILURES 


VARIABLE DEFINITIONS 


ALPHA:WEIBULL SHAPE PARAMETER 

AC:MATRIX USED TO STORE EACH AIRCRAFTS PREVIOUS FLIGHT HOURS 
WITHOUT FAILURE AND UNCONDITIONAL SURVIVAL FUNCTION 

AS:AIRCRAFT SELECTED FOR FLIGHT EVOLUTION 

CPF :CONDITIONAL PROBABILITY OF FAILURE 

DETMESDEPLOYMENT TIME 

٦.٣۱٣٦٣۴٣۸ USED TO STORE SIMULATED BALLURES AT SPECIFIED TIME 

INTERVALS FOR EACH DEPLOYMENT SIMULATION 

۶٠۳۷۶ ۶ ]ا‎ ٤5 

GFC:GENERATED FAILURE CRITERIA 

HM HISTORICAL MEAN 

HSD:HISTORICAL STANDARD DEVIATION 

LHAT:MAXIMUM LIKELIHOOD ESTIMATE FOR WEIBULL SCALE PARAMETER 

NAC:NUMBER OF AIRCRAFT USED IN SIMULATION 

NI:NUMBER OF TIME INTERVALS TO BE USED IN. SIMULATION (500 ds INT) 

HSENUHMBER- OB SIMULATIONS TO BE RUN 

OECD SEED USED TO DETERMINE AIRCRAFT CHOSEN FOR FLIGHT 

ETD SEED USED TO DETERMINE FLIGHT DURATION 

O EEDS SEED USED TO HDETERMINE IF A COMPONENT FAILURE WAS OBSERVED 

5626 ۰:06 GATE 

LIE ACCUMULATED FLIGHT TIME 

UCPS: UNCONDITIONAL PRORAILITY OF FAILURE 


DIMENSION AC(70,2),FAIL(20,100) ,EST(20,4) 


INTEGER AS 
REAL * 4 LHAT 


INITIALIZE 
DATA LHAT/. 0149/ALPHA/ . 8079 / 
DATA NI/16/NS/100/NAC/48/HM/2.18/HSD/ .9897/ 
DATA NSEED1/749684/NSEED2/4683957 /NSEED3/5872654/ 


COMMENCE SIMULATION(I) 
DO 1 5 
NF=0 
TIME=0. 
TGATE=500. 
DO 5 K=1,NAC 
AC p 
eli 251: 
CONTINUE 


SIMULATE AND RECORD COMPONENT FAILURES USING 500 HR INTERVALS 

DO 10 J=1,NI 
LRND AND LNORM ARE SUBROUTINES IN A RANDOM NUMBER GENERATOR PACKAGE 
CALLED LLRANDOMII: A NON-IMSL PROGRAM AVAILABLE AT THE NAVAL 
POSTGRADUATE SCHOOL. 


GENERATE UNIFORM (0,1) RANDOM VARIABLE TO DETERMINE WHICH AIRCRAFT 
WILL BE FLOWN 


CALL LRND(NSEED1,U,1,2,0) 
AS=(U*(NAC-1))+1.5 
GENERATE NORMAL (0,1) RANDOM VARIABLE TO COMPUTE FLIGHT DURATION 
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ېر DE DE ŽE RR R‏ پر پر پر پر پر پر ېز پر obo ob NE DE DE‏ پد پد 


Lt b خر پر‎ 


CALL LNORM(NSEED2,RN,1,2,0) 


FT=(RN*HSD)+HM 
TE e FT=.5 
EE ENEE 


COMPUTE THE UNCONDITIONAL PROBABILITY OF SURVIVAL FOR S+FT HRS, 
WHERE S EQUALS THE TOTAL FLIGHT HOURS FLOWN SINCE THE LAST 
COMPONENT FAILURE, AND FT EQUALS THE FLIGHT TIME OF THE SORTIE 
BEING FLOWN. 
UCPS=EXP(-1*((LHAT*(AC(AS,1)+FT))**ALPHA)) 
CPF=1-UCPS/AC(AS,2) 


= GENERATE UNIFORM (0,1) RANDOM VARIABLE TO DETERMINE IF COMPONENT 
x FAILURE OCCURS 


CALL LRND(NSEED3,GFC,1,2,0) 
* COMPONENT FAILURE p 
IF (GFC.LT.CPF) THEN 
NF=NF+1 
کت‎ 
AGUA pb. 
GO TO 15 
END IF 
k COMPONENT SURVIVES 


AC(AS, L)SAC(AS, 1) +FT . : 
AC(AS,2)<UCPS 


* UPDATE TIME 
15 TIME=TIME+FT 
IF (TIME.GE.TGATE) THEN 
FAIL(J,1)=NF 
TGATE=TGATE+500. 
COMO 10 
END IF 
* NEXT AIRCRAFT SORTIE TO BE SIMULATED 
GO TO l1 
* NEXT TIME INTERVAL 
10 CONTINUE 
NEXT DEPLOYMENT 
CONTINUE 
CALCULATE MEAN AND STANDARD DEVIATION FOR THE NUMBER OF FAILURES 
OBSERVED AT EACH TIME INTERVAL (MEAN=EST(J,1);STANDARD 
DEVIATION=EST(J,2)) 
DO 20 J=1,NI 
DO 25 I=1,NS 
E ES ک0‎ 
EST(J,2)=EST(J,2)+FAIL(J,I)*FAIL(J,I) 
25 CONTINUE 
20 CONTINUE 


DO 30 J=1,NI 


+ 4 x 4 


+ 


FA E 


100 


(I ESTO] Z2EST(J,1)))7 (NS=-1))%%*,5 
CONTINUE 


95 PERCENT Se INTERVAL (LOWER LIMIT=EST(J,3); 
UPPER LIMIT=EST(J,4) 


DO 35 J=1 
HINT=1. "O6*EST(J, 2) 
O, 3}= A 
EST(J,4)=EST(J,1)+HINT 
CONTINUE 


OUTPUT 


TIME=50 
WRITE 6 .1000( 
DO 40 UE Na 
WRITE (6, ۵ DEME TESE I I) EST(J,2) EST I3) EST(J,4) 
TIME=TIME+50 
CONTINUE 


30 


35 


40 


* OUTPUT FORMAT 


*1000 FORMAT ('1 TIME MEAN STD DEV LOWER UPPER') 
111١١ FORMAT (F5.0,4(2X,F7.3)) ۱ 


STOP 
END 
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Figure Li Flight Hour Models: Model Set WUC:73411100. 
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Figure 1.2 Geometric Model: Model Set WUC:734IT100. 
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Figure [.3 Flight Hour Models: Validation Set WUC:734H 100. 
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Figure 1.4 Geometric Model: Validation Set WUC:734H100. 
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Flight Hour Models: Model Set WUC:J60X 1600. 


Figure [.5 
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Figure 1.6 Geometric Model: Model Set WUC:46X 1600. 
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Figure 1.7 Flight Hour Models: Validation Set WUC:46X1600. 
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Figure 1.8 Geometric Model: Validation Set WUC:46x1600. 
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Figure 1.9 Flight Hour Models: Model Set WUC:4622100. 
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Figure 1.10 Geometric Model: Model Set WUC:4622100. 
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Figure 1.11 Flight Hour Models: Validation Set WUC:4622100. 
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Figure 1.12 Geometric Model: Validation Set WUC:4622100. 
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Figure 1.13 Flight Hour Models: Model Set WUC:56X 
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Figure 1.14 Geometric Model: Model Set WUC:56X 2500. 
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Figure 1.15 Flight Hour Models: Validation Set WUC:36X2500: 
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Figure 1.16 Geometric Model: Validation Set WUC:56X2500. 
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Figure 1.17 Flight Hour Models: Model Set WUC:632Z100. 
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Figure 1.18 Geometric Model: Model Set WUC:632Z100. 
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Figure 1.19 Flight Hour Models: Validation Set WUC:632Z100. 
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Figure 1.20 Geometric Model: Validation Set WUC:632Z100. 
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- Figure 1.22 Geometric Model: Model Set WUC:74A1500. 
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Figure 1.23 Flight Hour Models: Validation Set WUC:74A1500. 
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Figure 1.24 Geometric Model: Validation Set WUC:74A 1500. 
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Figure 1.25 Geometric Model: Model Set WUC:74A5M00. 
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Figure 1.26 Flight Hour Models: Validation Set WUC:74A5MO00. 
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Figure 1.27 Geometric Model: Validation Set WUC:74A5M00, 
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Figure 1.28 Flight Hour Models: Model Set WUC:5772200. 
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Figure 1.29 Geometric Model: Model Set WUC:5772200. 
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.30 Flight Hour Models: Validation Set WUC:5 
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Figure 1.31 Geometric Model: Validation Set WUC:5772200. 
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Figure 1.32 Flight Hour Models: Model Set WUC:6918100. 
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Figure 1.33 Geometric Model: Model Set WUC:6918100. 
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Figure 1.34 Flight Hour Models: Validation Set WUC:6918100. 
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Figure 1.35 Geometric Model: Validation Set WUC:691$ 100. 
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Ae foes Me AA AS E 6 " 
لئ یں‎ Z 



















ni i 


موا وا 






















هز اما 

A ابا‎ > 

لل Vy wy‏ لک ALA‏ یش پور ہنی نس سنا دمه ری نی نپ 

LL یا‎ leia d ex o Ir » 9v lp. "viti DR جاک کر با په . 1 ا‎ 
Een لکا نیا‎ ye ron ke Z doli dli alla ki ki di! Sai d, kr رو‎ > 

EN? aaa ai EE LIA AA EMT TEE DENK D TEEN d‏ خښ وسم سم په وا نا 

duda VAJE AN y 4 E‏ من Léen LATI II DI dli‏ ماب ها رنه له 05ر و vi ٢‏ سر 

di A KA AA رید یش یہ سیر ایی‎ AAA ČU, NT سيلا‎ 

PO DAN OPA A lt, echt Aer dén BA A A A ena Ae ېښ هس‎ À r 

tt a de dd de de AL A a li tea ai ULE e E Mech T 

Sse‏ هس نشی i a Deda A ADA RIA IE A DA Te ti EE farmi‏ سی8ز یلست لی 

A L. k. R 3 è‏ نیا زیر یا رر L Liia U e LIESTE T‏ له نت یج پر تی سے 

: پر تاور رھش رس تی یلال سلاد یس يل a JA RE Bel‏ ذف ون بی ہہ انت oA‏ 

dle di. ih ote Seg lil بک‎ ar jd کی کت یپ زیر یں ہر نپراز زار‎ ee ae H LA 1 TES a 5 s 

PREVI TUTE RARI PETE TTI A TT AMALIA L Pri y rd D eh شس‎ x. v. 

a O + VY AY toa da a DELA CA JA ha be T2 14 4 eu py یرہ‎ Rosly kA) 01 geg Vo A . : 

D ^‏ 2 ۲ ` رت پر لک GUANO Y us e RA‏ لهو یا did ida‏ هک نه هی اط بپ نف 

A II NA II PIAR y RA کہ و ارز لی‎ D We NEE SER s Te We AN? : مت‎ ٤ a 

ei dr TETI هن اه نا‎ Gr Se ett oe «v Pue RI کر جج‎ ` 7 ' 

dd ide ad ia o A A La سینا‎ LL HC Ta | "T í e ۲ 

V4 S8 (79 Cp lrf ro nd ax prot FA nan P CADE COM P T هد دە هه‎ an Š « کی‎ ` : 

e هغمه مهد رس لي‎ OD rale e qh LP Léa ln An Q چس‎ 

y ay Ji Li 7 ENE T m a d . d 
Ebro (sto Ja Los risit Bed بب‎ Jc ld LE 7 رہ اا‎ H e pb Y t s: Li V è ñ % 
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